2202.05892v2 [astro-ph.SR] 7 Aug 2022 


arXiv 


DRAFT VERSION AUGUST 9, 2022 
Typeset using IATEX twocolumn style in AASTeX63 


POSYDON: A General-Purpose Population Synthesis Code with Detailed Binary-Evolution Simulations 


Tassos FRAGOS,! JEFF J. ANDREWS,”? SIMONE S. BAVERA,! CHRISTOPHER P.L. BERRY,?* SCOTT COUGHLIN,” 
AARON DOTTER,” PRABIN GIRI,” VICKY KALOGERA,? AGGELOS KATSAGGELOS,” " KONSTANTINOS KOVLAKAS,! 
SHAMAL LALVANI,’ DEVINA MisRA,! PHILIPP M. SRIVASTAVA,” ” YING Qin,® KYLE A. Rocna,”® JAIME ROMAN-GaARZA,® 
JUAN GABRIEL SERRA,” PETTER STAHLE, MENG SuN,? Xu TENG,’ GOCE TRAJCEVSKI,? NAM Har TRAN,!! ZEPEI XING,! 
EMMANOUIL ZAPARTAS,! AND MICHAEL ZEvIN?: 8 


1 Département d’Astronomie, Université de Genève, Chemin Pegasi 51, CH-1290 Versoix, Switzerland 


? Center for Interdisciplinary Exploration and Research in Astrophysics (CIERA), Northwestern University, 2145 Sheridan Road, 
Evanston, IL 60208, USA 


3 Department of Physics, University of Florida, 2001 Museum Rd, Gainesville, FL 32611, USA 
4 Institute for Gravitational Research, University of Glasgow, Kelvin Building, University Avenue, Glasgow, G12 8QQ, Scotland 
5 Department of Electrical and Computer Engineering, Iowa State University, 2520 Osborn Dr, Ames, IA 50011, USA 
6 Department of Physics and Astronomy, Northwestern University, 2145 Sheridan Road, Evanston, IL 60208, USA 
" Electrical and Computer Engineering, Northwestern University, 2145 Sheridan Road, Evanston, IL 60208, USA 
8 Department of Physics, Anhui Normal University, Wuhu, Anhui 241000, China 
9 Universidad de Monterrey, Ave. Morones Prieto 4500 Pte., C.P. 66283, San Pedro Garza García, Nuevo León, México. 
10 Département d'Informatique, Université de Genève, Route de Drize 7, CH-1227 Carouge, Switzerland 
11 DARK, Niels Bohr Institute, University of Copenhagen, Jagtvej 128, 2200 Copenhagen, Denmark 
1? Kavli Institute for Cosmological Physics, The University of Chicago, 5640 South Ellis Avenue, Chicago, Illinois 60637, USA 
13 Enrico Fermi Institute, The University of Chicago, 933 East 56th Street, Chicago, Illinois 60637, USA 


Submitted to AAS Journals 


ABSTRACT 


Most massive stars are members of a binary or a higher-order stellar systems, where the presence 
of a binary companion can decisively alter their evolution via binary interactions. Interacting binaries 
are also important astrophysical laboratories for the study of compact objects. Binary population 
synthesis studies have been used extensively over the last two decades to interpret observations of 
compact-object binaries and to decipher the physical processes that lead to their formation. Here, we 
present POSYDON, a novel, publicly available, binary population synthesis code that incorporates full 
stellar-structure and binary-evolution modeling, using the MESA code, throughout the whole evolution 
of the binaries. The use of POSYDON enables the self-consistent treatment of physical processes in 
stellar and binary evolution, including: realistic mass-transfer calculations and assessment of stability, 
internal angular-momentum transport and tides, stellar core sizes, mass-transfer rates and orbital 
periods. This paper describes the detailed methodology and implementation of POSYDON, including the 
assumed physics of stellar- and binary-evolution, the extensive grids of detailed single- and binary-star 
models, the post-processing, classification and interpolation methods we developed for use with the 
grids, and the treatment of evolutionary phases that are not based on pre-calculated grids. The first 
version of POSYDON targets binaries with massive primary stars (potential progenitors of neutron stars 
or black holes) at solar metallicity. 


1. INTRODUCTION 


Throughout their lives, stars affect their surroundings 
via the immense energy radiated across the electromag- 


netic spectrum (e.g., Conroy 2013; Eldridge & Stanway 
2022) and the nuclear-processed material emitted as a 
stellar wind (e.g., Kudritzki & Puls 2000; Smith 2014). 
The deaths of massive (>, 8 Mọ) stars, even more than 
their lives, transform their environments as their cores 
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run out of nuclear fuel and collapse to form neutron 
stars (NSs) and black holes (BHs). The formation of 
these compact objects (COs) is often accompanied by a 
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supernova (SN) or a y-ray burst that release more en- 
ergy in 10s than our Sun in 10!° yr (e.g., Janka 2012; 
Burrows & Vartanyan 2021; Woosley & Bloom 2006). 
These explosive events enrich their environments with 
heavier elements while also regulating any ongoing star- 
formation (e.g., Nomoto et al. 2013; Hopkins et al. 2011, 
2012). 

It is now established that most massive stars are mem- 
bers of a binary or a higher-order stellar system (Sana 
et al. 2013; Moe & Di Stefano 2017). More often than 
not, the presence of a binary companion decisively al- 
ters the evolution and final fate of both binary compo- 
nents via binary-interaction processes such as tidal dissi- 
pation, mass-transfer phases, and stellar mergers (Sana 
et al. 2012; De Marco & Izzard 2017). Furthermore, in- 
teracting binaries are arguably some of the most impor- 
tant astrophysical laboratories available for the study 
of COs. Accretion of matter from a binary companion 
gives rise to X-ray emission, bringing the system to the 
X-ray binary (XRB) phase (Bhattacharya & van den 
Heuvel 1991; Podsiadlowski et al. 1992), while gravi- 
tational waves (GWs) enable us to witness the last mo- 
ments of the lives of coalescing binary COs (Abbott et al. 
2016, 2021). 

Over the last two decades, multi-wavelength surveys 
of the Milky Way and its neighborhood, as well as nu- 
merous nearby and more distant galaxies, have amassed 
large datasets of binary stellar systems. These datasets 
range from targeted, high-resolution observations of 
galaxies in the local Universe (e.g., the PHAT and LE- 
GUS surveys; Dalcanton et al. 2012; Calzetti et al. 2015), 
to serendipitous (e.g., the Chandra Source Catalogue 
2.0; Evans et al. 2019), all-sky (e.g., Gata; Gaia Col- 
laboration et al. 2018) and transient surveys (e.g., Pan- 
STARRS, Zwicky Transient Facility; Kaiser et al. 2002; 
Graham et al. 2019). In addition to electromagnetic 
surveys, there is also the global GW observatory net- 
work of LIGO (Aasi et al. 2015), Virgo (Acernese et al. 
2015) and KAGRA (Akutsu et al. 2019), which have 
detected nearly a hundred binary CO mergers (Abbott 
et al. 2021). Combined, these surveys are revolution- 
izing our view of binary stellar systems, including CO 
binaries, and their environments. 

Aspects of the astrophysics of all these different types 
of stellar binaries can be obtained from observations and 
modeling of present-day properties of individual, well- 
studied systems. However, more comprehensive insight 
requires understanding the statistical properties of their 
entire populations. For these studies, binary population 
synthesis (BPS) modeling is often employed. BPS mod- 
eling first generates initial binary populations, whose 
properties are randomly sampled from probability distri- 


butions that can be observationally constrained. Then, 
this initial population is evolved with a computationally 
efficient simulation tool using our best understanding of 
the physics dictating binary star interactions, to pro- 
duce observable properties of the target population. If 
the number of binaries evolved is large enough to provide 
a statistically significant description of a population of 
interest, then BPS can provide valuable insights about 
the expected rate and distribution of the target popu- 
lation's properties, the different evolutionary pathways 
that lead to formation of these systems, and the effect 
that different physical processes have on their evolution. 

Over the last two decades, many general purpose 
BPS codes have been developed, e.g., binary_c (Izzard 
et al. 2004, 2006, 2009), BPASS (Eldridge et al. 2017), 
the Brussels code (Vanbeveren et al. 1998a,b), BSE 
(Hurley et al. 2002), ComBinE (Kruckow et al. 2018), 
COMPAS (Stevenson et al. 2017; Riley et al. 2022), COSMIC 
(Breivik et al. 2020), MOBSE (Giacobbo et al. 2018), the 
Scenario Machine (Lipunov et al. 1996, 2009), SEVN 
(Spera et al. 2015), SeBa (Portegies Zwart & Verbunt 
1996; Toonen et al. 2012a), StarTrack (Belczynski et al. 
2002, 2008), and TRES (Toonen et al. 2016). These have 
been used in studies of a wide variety of binary popu- 
lations. A hard requirement for BPS is computational 
efficiency, as for most studies one would need to model 
the evolution of many millions of binaries in a reasonable 
computational time. 

BPS codes stand in stark contrast to detailed stellar- 
structure and binary-evolution codes, e.g., BEC (Heger 
et al. 2000; Heger & Langer 2000), BINSTAR (Siess et al. 
2013), the Cambridge STARS code (Eggleton 1971; Pols 
et al. 1995; Eldridge & Tout 2004; Stancliffe & Eldridge 
2009), MESA (Paxton et al. 2015), and the TWIN code 
(Nelson & Eggleton 2001; Eggleton & Kiseleva-Eggleton 
2002), which self-consistently solve the stellar struc- 
ture equations of a binary's component stars along with 
the orbital evolution. Many studies have used detailed 
binary-evolution calculations (e.g., Nelson & Eggleton 
2001; Podsiadlowski et al. 2002; de Mink et al. 2007; 
Marchant et al. 2017; Qin et al. 2018, 2019; Langer et al. 
2020; Misra et al. 2020; Laplace et al. 2020, 2021) to 
generate grids of models, varying the masses of the two 
stars and the binary's orbital period. However, in all 
those cases, the grids of detailed binary tracks either 
cover a limited part of the initial parameter space, or 
focus on a specific evolutionary phase. This limitation 
is principally caused by the computational demands of 
detailed grids; each simulation typically requires ~ 10- 
100 CPU hours for the modeling of a single system (e.g., 
Paxton et al. 2019). 
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As a result of this computational expense, a common 
thread among the vast majority of current BPS codes 
is that they approximate each star’s evolution, employ- 
ing either fitting formulae (e.g., SSE; Hurley et al. 2000) 
or look-up tables (e.g., COMBINE; Kruckow et al. 2018) 
for the properties of single stars, based on grids of pre- 
calculated detailed, single-star models. Then, the ef- 
fects of binary interactions (e.g., Roche-lobe overflow 
or tides) are modeled using approximate prescriptions 
and parametrizations. This modeling approach is often 
called rapid or parametric BPS; throughout the remain- 
der of this work, we choose to use the term parametric 
BPS (pBPS) modeling, to make a distinction between 
computational efficiency and modeling accuracy. A no- 
table exception among BPS codes is BPASS (Eldridge 
et al. 2017), which uses extensive grids of detailed binary 
evolution models computed with a custom version of the 
Cambridge STARS binary evolution code (Stancliffe & 
Eldridge 2009). In the grids of binary-star models em- 
ployed in BPASS, both the primary and the secondary 
stars are followed in detail, but only one at a time (for 
computational-cost reasons). During the primary’s evo- 
lution, the properties of the secondary star are approxi- 
mated by formulae based on single-star models (Hurley 
et al. 2000). Subsequently, once the modeling of the 
primary’s evolution is completed, the secondary star’s 
evolution is re-computed, accounting for mass-transfer 
and rejuvenation effects. 

Studies using pBPS techniques have allowed us to 
make advances in our understanding of the formation 
pathways leading to different types of binary systems, 
and to interpret observations of binary populations. 
Such examples are studies on white-dwarf binaries (e.g., 
Nelemans et al. 2001; Ruiter et al. 2009, 2010; Too- 
nen et al. 2012b; Breivik et al. 2018; Korol et al. 2020), 
XRBs (e.g., Van Bever & Vanbeveren 2000; Belczynski 
et al. 2004; Fragos et al. 2008, 2013b,a; Luo et al. 2012; 
Tzanavaris et al. 2013; Tremmel et al. 2013; Zuo et al. 
2014; Zuo & Li 2014; van Haaften et al. 2015; Artale 
et al. 2019; Wiktorowicz et al. 2019; Shao & Li 2020), the 
Galactic population of double NSs (DNSs; e.g., Ostowski 
et al. 2011; Andrews et al. 2015; Chruslinska et al. 2017; 
Vigna-Gómez et al. 2018; Chattopadhyay et al. 2020), 
GW sources observable by ground-based observatories 
(e.g., Dominik et al. 2012, 2013, 2015; Mennekens & 
Vanbeveren 2014, 2016; Belczynski et al. 2016; Klencki 
et al. 2018; Giacobbo et al. 2018; Mapelli & Giacobbo 
2018; Neijssel et al. 2019; Spera et al. 2019; Kinugawa 
et al. 2020; Breivik et al. 2020; Zevin et al. 2020; Broek- 
gaarden et al. 2021) and SNe in binary systems (e.g., 
De Donder & Vanbeveren 2003; Vanbeveren et al. 2013; 
Claeys et al. 2014; Zapartas et al. 2017, 2019, 2021). 


However, the implicit assumption in pPBS codes that 
the binary components have properties identical to sin- 
gle stars of the same mass in thermal equilibrium (e.g., 
abundance profiles, core sizes, mass-radius relations and 
response to mass-loss), as well as the lack of information 
about the star's internal structure at different critical 
evolutionary phases (e.g., the onset of a dynamically 
unstable mass-transfer, the end of stable and unstable 
mass-transfer phases or the core-collapse), may intro- 
duce systematic uncertainties and inaccuracies. Cur- 
rent pBPS codes therefore rely on approximate prescrip- 
tions for modeling binary interactions and difficult-to- 
calibrate additional model parameters. These compli- 
cations could be avoided by instead employing detailed 
stellar-structure and binary-evolution simulations (here- 
after detailed models). Focusing on aspects that are rel- 
evant to the formation of CO binaries, detailed models 
(i) allow for a self-consistent estimation of the mass- 
transfer rate, especially during thermal-timescale mass- 
transfer phases, and therefore an accurate assessment 
of mass-transfer stability; (ii) allow for a more accurate 
description of the type and properties of the formed CO 
as well as any potential associated transient events since 
the internal structure of pre-core-collapse stars is known, 
(iii) account for the transport of angular momentum 
between and within the binary components, including 
its back-reaction on the structure and evolution of each 
star (e.g., rotational mixing), and (iv) allow for the self- 
consistent modeling of the end of a mass transfer phase 
(e.g., accounting for a potential partial stripping of the 
envelope). 

In this work, we build upon the combined experience 
gained from the large body of BPS studies to date, to 
create POSYDON (POpulation SYnthesis with Detailed 
binary-evolution simulatiONs), a general-purpose code 
that can generate entire populations of binaries, under- 
pinned by detailed, self-consistent models of stellar bi- 
naries.! With POSYDON we aim to address many of the 
caveats of pBPS codes, while at the same time maintain- 
ing much of their flexibility. In its first release (v1.0), 
POSYDON is limited to stars of solar metallicity, and bi- 
naries where the primary star is massive enough to form 
a BH or a NS. Future releases, which are already in de- 
velopment, will lift these limitations. In Section 2, we 
introduce POSYDON and the approach it takes to model- 
ing binary populations. In Sections 3 and 4 we provide 
the physics adopted for our detailed models of single 
and binary stars, respectively. We describe the pre- 


1 POSYDON will become is publicly available at https://posydon.org 
blieati Cg 
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calculated grids of single and binary stellar evolution 
models in Section 5, the way they are post-processed in 
Section 6, and our classification and interpolation meth- 
ods for their optimal use in Section 7. In Section 8 
we detail our treatment of evolutionary phases which 
are not based on pre-calculated grids, such as the core- 
collapse and the common-envelope (CE) phase, while in 
Section 9 we describe how all the aforementioned pieces 
come together to model the entire evolution of a binary. 
In Section 10 we outline our assumptions and methods 
in modeling populations of binary systems and present 
some example results. We conclude in Section 11, where 
we present an outlook of future development directions 
of the POSYDON code. The definitions of all the symbols 
used throughout this paper can be found in Table 1. 


Table 1. List of variables used throughout the paper. 


Name Description First appears 
a Orbital separation 4.1 
ai Orbital separation before orbital 8.3.5 
kick 

af Orbital separation after orbital 8.3.5 
kick 

Üpre,CE Orbital separation pre common 8.2 
envelope 

Opost,CE Orbital separation post common 8.2 
envelope 

a Rate of change of orbital 8.1.2 
separation 

Awind Rate of change of orbital separa- 8.1.2 
tion due to wind mass loss 

Atides Rate of change of orbital separa- 8.1.2 
tion due to tides 

üGR Rate of change of orbital sepa- 8.1.2 
ration due to gravitational-wave 
radiation 

aspin Non-dimensional spin 4.2.3 

c Speed of light 4.2.2 

Dconv.reg. Depth of a convective region TA 

[2 Orbital eccentricity 8.1.2 

ef Orbital eccentricity after orbital 8.3.5 
kick 

é Rate of change of orbital 8.1.2 
eccentricity 

Ctides Rate of change of orbital eccentric- 8.1.2 
ity due to tides 

ÉGR Rate of change of orbital eccen- 8.1.2 
tricity due to gravitational-wave 
radiation 

E Eccentric anomaly 8.3.5 

E» Second-order tidal torque 4.1 
coefficient 


'Table 1 continued 


Table 1 (continued) 


Name Description First appears 
fiony Dimensionless factor accounting 4.1 
for slow convective shells that can- 
not contribute to the tidal viscosity 
within an orbital timescale 
ft» Fallback mass fraction 8.3.2 
foy Convective exponential overshoot- 3.2.3 
ing parameter 
g Local gravitational acceleration 3.2.1 
G Gravitational constant 3.2.2 
I Moment of inertia 4.1 
" Moment of inertia rate of change 8.1.2 
j Specific angular momentum 5.7 
jisco Specific angular momentum of the 8.3.4 
ISCO 
J Angular momentum 5.7 
Jshell Stellar shell’s angular momentum 8.3.4 
Jairect Stellar shell’s angular momentum 8.3.4 
of directly collapsing material 
Jaisk Stellar shell’s angular momentum 8.3.4 
of disk forming material 
JBu Black hole angular momentum 8.3.4 
k Apsidal motion constant 4.1 
L Star's luminosity 3.2.2 
Lgaa Eddington luminosity 3.2.2 
Le Second lagrange point 8.2 
Mi Mass of the initially more massive 5.5 
star 
Mo Mass of the initially less massive 5.5 
star 
Macc Mass of the accretor 4.2.2 
Maon Mass of the donor 8.2 
Mco Mass of the compact object 5.6 
Moj;O-core Mass of the C/O core 7.4 
Mconv.reg. Mass of a convective region 4.1 
Mcomp Mass of the binary companion star 4.1 
Maisk Mass of accretion disk 8.3.4 
Menv Mass of the stellar envelope 7.4 
Mgrav Remnant’s gravitational mass 8:9.3 
Myte-core Mass of the He core 7.4 
Miembar Remnant's baryonic mass 7.4 
MNS” Maximum neutron star mass 8.3.3 
Mgaa Mass-accretion rate corresponding 4.2.2 
to the Eddington limit 
My Wind mass-loss rate 92.2.2 
Mi ot Binary stellar mass before core 8.3.5 
collapse 
Mí. Binary stellar mass after core 8.3.5 
collapse 
Meshell Stellar shell’s mass 8.3.4 
P Pressure 3.2.1 


Table 1 continued 
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Table 1 (continued) 


Name Description First appears 

Fi Orbital period 5.5 

q Binary mass ratio 4.1 

R Stellar radius 3.2.2 

Race Radius of the accretor 4.2.2 

FHi,conv.reg. Radial coordinate of a convective 4.1 
region’s bottom boundary 

Reore Radius of the stellar core 7.4 

Ro/o-core Radius of the C/O core 7.4 

Reonv Radius of the convective core 4.1 

Reonv.reg. Radial coordinate of a convective 7.4 
region’s center 

Raon Radius of the donor 8.2 

Rye—core Radius of the He core 7.4 

Ri, Roche lobe radius 4.2.2 

Ry ace Roche lobe radius of the accretor 4.2.2 

Rt,conv.reg. Radial coordinate of a convective 4.1 
region’s top boundary 

r Stellar shell’s radius 8.3.4 

p Instantaneous orbital separation 8.3.5 
before orbital kick 

Tef Effective temperature 3.2.1 

T Timescale for orbital changes due 4.1 
to tides 

Uk Magnitude of velocity kick 8.3.5 

Ur Orbital velocity of the collapsing 8.3.5 
star 

X Hydrogen mass function 4.2.2 

A center Center hydrogen mass function 8.1.1 

Xsurf Surface hydrogen mass function S.I 

ba Helium mass fraction 3.1 

Voenter Center helium mass fraction 8&1. 

Yourk Surface helium mass fraction 7.5 

Z Metallicity (mass fraction of ele- 3.1 
ments heavier than ^He) 

«GE Fraction of the orbital energy that 8.2 
contributes to the unbinding of the 
CE 

QMLT Convective mixing length 3.2.3 
parameter 

Oth Thermohaline mixing parameter 3.2.8 

n Dimensionless factor denoting the 4.2.2 
radiative efficiency of the accretion 
process 

0 Stellar profile’s polar angle 8.3.4 

Odisk Stellar profile's polar angle of disk 8.3.4 
formation 

K Opacity 3.2.1 

ACE Parametrization of the CE’s bind- 8.2 
ing energy 

T Optical depth 3.2.1 

Tsync 'Tidal synchronization timescale 4.1 


'Table 1 continued 


Table 1 (continued) 


Name Description First appears 
Tony Convective timescale 4.1 
Tmb Magnetic braking torque 8.1.2 
V Reduced mass 8.1.2 
c Maxwellian distribution dispersion 8.3.5 
OCCSN Maxwellian distribution dispersion — 8.3.5 
for CCSN kicks 

CECSN Maxwellian distribution dispersion 8.3.5 
for ECSN kicks 

ab Binary orbital inclination with re- 8.3.5 
spect to before the kick 

Us Surface angular velocity 3.2.2 

Us, crit Critical surface angular velocity 

Qorb Orbital angular velocity 8.1.2 

OQshell Stellar shell’s angular velocity 8.3.4 

Q Stellar angular velocity 8.1.2 

Q Stellar angular velocity rate of 8.1.2 
change 

ud Stellar angular velocity rate of 8.1.2 
change due to winds 

Cisestta Stellar angular velocity rate of 8.1.2 
change due to changes of star's mo- 
ment of inertia 

Prides Stellar angular velocity rate of 8.1.2 
change due to tides 

mb Stellar angular velocity rate of 8.1.2 
change due to magnetic breaking 3.2.2 


2. OVERVIEW OF THE STRUCTURE OF POSYDON 


At its core, a BPS code requires two elements: a 
method to generate random binaries at zero-age main 
sequence (ZAMS) and a mechanism to evolve those bi- 
naries. The former is described in Section 10.1. Regard- 
ing the evolution of each binary, our primary goal with 
POSYDON is to self-consistently evolve the internal struc- 
tures of the two stars comprising a binary along with 
the binary’s orbit. To achieve this goal, we have opted 
to employ the stellar-structure code Modules for Exper- 
iments in Stellar Astrophysics (MESA; Paxton et al. 2011, 
2013, 2015, 2018, 2019). However, calculating the evo- 
lutionary track of a binary using MESA can take in excess 
of 100 CPU hours, a six-orders-of-magnitude increase in 
computational cost compared to a typical pBPS code 
such as COSMIC (Breivik et al. 2020). This means that, 
even with modern computing resources, we cannot fea- 
sibly run more than ~ 10°—10® binaries, far too few 
to accurately model a Milky Way sized population with 
~ 10! stellar binaries. As a further complication, codes 
like MESA cannot run individual binaries from start to 
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finish; key physics, including CE phases and supernovae 
require the code to be stopped and restarted.” 

POSYDON solves the problems associated with stellar- 
structure codes by using extensive, pre-calculated grids 
of single and binary stellar-evolution models, covering 
the parameter space relevant for the formation of high- 
mass binary stars, with separate grids being calculated 
for each phase of binary evolution. In v1.0 our grids 
contain a combined total of nearly 120,000 separate 
detailed binary simulations. To compute these grids, 
POSYDON has an infrastructure specifically designed for 
high-performance computing environments that stream- 
lines the process of producing large grids with consistent 
physics inputs. Using this infrastructure, we have gen- 
erated five separate grids of single and binary stellar 
evolution models. We computed three grids of inter- 
acting binary stars initially composed of two hydrogen 
(H)-rich ZAMS stars (Section 5.5); a CO and a H-rich 
star at the onset of Roche-lobe overflow (RLO; Section 
5.6), and a helium (He)-rich ZAMS star with a CO com- 
panion (Section 5.7). We further computed two grids of 
single H-rich and He-rich stars (Sections 5.3 and 5.4, re- 
spectively), which we use for the modeling of detached, 
non-interacting binaries. These five grids are then post- 
processed, so that their data size is reduced (Section 6). 
We additionally apply classification and interpolation 
algorithms on the outputs of these extensive grids (Sec- 
tion 7), allowing us to effectively interpolate between 
MESA simulations to estimate the evolution of any ar- 
bitrary binary within some bounded region of the pa- 
rameter space. As a simpler alternative, we also provide 
functionality to evolve individual binaries using nearest- 
neighbor matching, and in that case no classification or 
interpolation methods are required. 

The second major component of POSYDON is the code 
infrastructure to follow the entire evolution of a binary 
from start to end. To achieve this, we combined the 
aforementioned grids (and classification and interpola- 
tion methods) with physics dictating a binary’s evolu- 
tion through key phases, including core collapse and CO 
formation (Section 8.3) and CE (Section 8.2). These 
latter phases are not modeled based on pre-calculated 
grids, but rather with on-the-fly calculations. Simi- 
larly, the evolution of detached, potentially eccentric, 
post-core-collapse binaries is also modeled with on-the- 
fly calculations, where we use the single-star grids cou- 
pled to binary evolution routines, i.e. orbital evolu- 


tion due to tides, stellar winds, magnetic breaking and 
gravitational-wave emission (Section 8.1). 

In the POSYDON approach, each separate evolution- 
ary phase (step) has its own dedicated function which 
determines the binary's state resulting from that step, 
the quantitative values characterizing that binary (e.g., 
masses of the two stars), and the event describing how 
that state ended (e.g., onset of RLO). Once a step is 
completed, the POSYDON framework uses the resulting 
binary state and event as well as each component stars' 
states to determine an individual binary's next evolu- 
tionary step. The process is repeated until a binary's 
evolution is complete, resulting in a disrupted binary, a 
binary merger, or a double CO. At this point, the next 
binary is run. The modular nature of POSYDON allows 
a user to also provide their own prescriptions to model 
each phase of evolution, or even their own breakdown 
of the binary-evolution tree. As default in POSYDON, we 
provide a complete set of evolutionary steps which we 
visually summarize in Figure 1 and present in detail in 
the following sections. 


3. ADOPTED STELLAR PHYSICS 


All stellar-evolution models described in this paper 
were computed with the state-of-the-art, open source 
stellar-structure and evolution code MESA (Paxton et al. 
2011, 2013, 2015, 2018, 2019) revision 11701 together 
with the 20190503 version of the MESA software devel- 
opment kit (SDK; Townsend 2020). MESA solves the 
one-dimensional stellar-structure and composition equa- 
tions. Mixing and burning processes are solved simulta- 
neously; mixing is treated as a diffusive process. Discus- 
sion of specific elements of stellar physics are described 
in the following subsections, which are split into micro- 
physical and macrophysical processes. We implement 
all physics that are not readily available in MESA us- 
ing the functionality provided by run star extras and 
run binary.extras. 


3.1. Microphysics 


We adopt the Asplund et al. (2009) protosolar abun- 
dances as our initial composition, with Z — 0.0142 and 
Y — 0.2703. The equation of state is the standard MESA 
amalgamation of the SCVH (Saumon et al. 1995), OPAL 
(Rogers & Nayfonov 2002), HELM (Timmes & Swesty 
2000) and PC (Potekhin & Chabrier 2010) equations 
of state (Paxton et al. 2019). Radiative opacities are 


? In principle these phases could be run within a stellar-structure 
code like MESA; for instance, more updated versions of MESA than 
the one we use can handle the evolution of a binary through a 
CE (Marchant et al. 2021). 


3 We made one minor bug fix in the MESA source code, which in- 
volves replacing the mass of the proton with the atomic mass 
unit where it appears in the code that evaluates the Potekhin & 
Chabrier (2010) equation of state (E. Bauer, 2020, private com- 
munication). This change is included in later MESA releases. 
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Figure 1. The structure of POSYDON (v1.0) for modeling the evolution of a binary star. Rectangles represent the initial and 
possible final outcomes of the evolution, and black circles represent events in the evolution of a binary. Colored lines showcase 
the different evolutionary steps that POSYDON follows. Evolutionary steps that are based on pre-calculated grids of detailed 
binary-evolution tracks are designated with solid lines, while those that are based on computations performed on-the-fly for 


each modeled binary are shown with dashed lines. 


taken from Ferguson et al. (2005) and Iglesias & Rogers 
(1996) for the Asplund et al. (2009) mixture, along with 
electron conduction opacities from Cassisi et al. (2007). 
Nuclear reaction rates are drawn from the JINA Rea- 
clib database (Cyburt et al. 2010). All models were 
computed using the approx21 nuclear reaction network 
that consists of 21 species: !H, ?He, *He, !?C, 14N, 160, 
20Ne, 7^ Mg, 285i, 325, 36 Ar, Os, Ti, 48Cr, 56Cr, 52Fe, 
54Fe, 56Fe and °°Ni, plus protons and neutrons (for the 
purpose of photo-disintegration). 


3.2. Macrophysics 
3.2.1. Surface Boundary Conditions 


The stellar surface boundary condition is satisfied by 
the simple-photosphere option, which sets the photo- 
sphere temperature using the Eddington (1926) Teg (T) 
relation, and the photosphere pressure via P = Tg/k 
with enhancement due to radiation pressure at the pho- 
tosphere (e.g., Paxton et al. 2011). 


3.2.2. Stellar Winds 


Stellar winds are a complex subject due to the varied 
physical mechanisms (known and unknown) that drive 
them and their dependence on the evolutionary state of 
their parent star. With this in mind, we have kept the 


wind prescription as simple as possible, while still cap- 
turing the key phenomenology of massive stellar evolu- 
tion, but avoiding fine-tuning to reproduce any single 
subset of observations. In POSYDON, changing the wind 
prescription would require the computation of new set 
of single- and binary-star model grids. 

For stars with initial masses above 8 Mo; we use the 
MESA Dutch scheme, which consists of de Jager et al. 
(1988) for Teg < 10,000K and Vink et al. (2001) for 
Tog > 11,000 K. In cases where Tig > 11,000 K and the 
surface 'H mass fraction is below 0.4, the Vink et al. 
(2001) wind is replaced with the Wolf-Rayet wind of 
Nugis & Lamers (2000). Between 10,000 K and 11,000 K 
there is a linear ramp (as a function of Teg) between the 
two wind prescriptions. We do not explicitly include any 
luminous blue variable (LBV) type winds; however, our 
stellar models at solar metallicity (e.g., Figure 3) do not 
enter the regime where LBV-type winds are typically 
applied in other studies (e.g., Belczynski et al. 2010). 

For stars with initial masses below 8 Mo, we again 
use the Dutch scheme for stars with Tig hotter than 
12,000 K. For stars with Tog less than 8,000 K, we use 
the Reimers (1975) wind with scaling factor nr = 0.1 
for stars on the first ascent of the giant branch, and the 
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Bloecker (1995) wind with scaling factor ng = 0.2 for 
stars in the thermally-pulsating phase. For the case of 
8,000K < Tog < 12,000K, we calculate the wind for 
both the hot and cool schemes, and linearly interpolate 
between the two. 

For mass-loss rates that have an explicit dependence 
on metallicity Z, we rescale wind mass loss based on 
the initial metallicity, not the current surface Z, as 
winds are driven predominantly by iron-group elements 
that remain almost constant throughout stellar evolu- 
tion (e.g., Vink & de Koter 2005). The primary motiva- 
tion for this approach is to avoid the dredge-up of carbon 
and oxygen to the surface layers in the later phases of 
evolution, which can cause surface Z to approach 1, from 
unduly influencing the mass-loss rate. The only excep- 
tion here is the wind prescription by Nugis & Lamers 
(2000) for Wolf-Rayet stars, which is specifically cali- 
brated to the total surface metal content, including car- 
bon and oxygen: in this case, we use the current, surface 
Z value of the He-rich star. 

We further boost stellar winds to limit a star’s rotation 
below its critical threshold (ws/ws,crit € 1), so that the 
sum of the centrifugal force and the photon pressure 
never exceeds gravity on the surface of the star. The 
impact of rotation on the mass loss rate is considered as 
indicated in (Heger & Langer 1998; Langer 1998), 


inoino) 


where Mw is the star’s wind mass-loss rate, and ws 
and Ws crit are the angular velocity and critical angu- 
lar velocity at the surface, respectively. The default 
value of the exponent € = 0.43 is taken from Langer 
(1998). The critical angular velocity is given the ex- 
pression w? i, = (1 — L/Leaa)GM/R*, where Lgaa is 
the Eddington luminosity and its expression is given in 
Eq. (8). This explicit boost to the wind is supplemented 
by an implicit numerical scheme implemented in MESA 
which ensures that the rotation of a star never exceeds 
its critical value. 


3.2.3. Convection, Rotation, and Mixing Processes 


Convective energy transport is modeled using mixing 
length theory (MLT; Bóhm-Vitense 1958) except in su- 
peradiabatic, radiation-dominated regions where we em- 
ploy the MLT++ modifications introduced in MESA (Pax- 
ton et al. 2013) that reduce the superadiabaticity in 
radiation-dominated convective regions, to improve nu- 
merical convergence. For the condition of convective 
neutrality we use the Ledoux criterion, and we use the 
convective premixing scheme as described by Paxton 
et al. (2019). We adopt a solar-calibrated mixing length 


parameter, au@iT = 1.93, based on results from the MIST 
project (Dotter et al., in preparation). 

Rotation is implemented in MESA as described in Pax- 
ton et al. (2013, 2019). Rotational mixing and angular- 
momentum transport follow the MIST project (Choi 
et al. 2016). It has been suggested that magnetic 
angular-momentum transport processes are main candi- 
dates for efficient coupling between the stellar core and 
its envelope during the post-MS (post-main sequence). 
Here, we adopt the Spruit-Tayler (ST) dynamo (Spruit 
2002) that can be produced by differential rotation in 
the radiative layers and amplify a seed magnetic field. 
Stellar models with ST dynamo can reproduce the flat 
profile of the Sun (Eggenberger et al. 2005) and observa- 
tions of the final spins of both white dwarfs and neutron 
stars (Heger et al. 2005; Suijs et al. 2008), but struggles 
to explain the slow rotation rates of cores in red giants 
(Eggenberger et al. 2012; Cantiello et al. 2014; Fuller 
et al. 2019). 

MESA treats mixing processes in the diffusive approx- 
imation with MLT providing the basic description. In 
addition to MLT convection, we consider thermohaline 
mixing with the parameter ay, = 17.5 (Paxton et al. 
2013, Eq. 14), also referred to as C4 (Charbonnel & Zahn 
2007, Eq. 4), corresponding to an aspect ratio ~ 1 of the 
instability fingers (Kippenhahn et al. 1980). Thermoha- 
line mixing is important during mass accretion from an 
evolved primary star onto an unevolved secondary star 
because the accreting material typically has a higher 
mean molecular weight than the material near the sur- 
face of the secondary star (e.g., Kippenhahn et al. 1980). 

Overshoot mixing is treated in the exponential de- 
cay formalism (Herwig 2000; Paxton et al. 2011). For 
the parameter foy describing the extent of the over- 
shoot mixing in this formalism, we adopt an initial- 
mass-dependent relation. For lower mass stars (initial 
masses less than 4 Mo) we adopt a value taken from 
the MIST project, fo, = 0.016, which is calibrated us- 
ing the Sun, as well as open clusters (Choi et al. 2016). 
In the high-mass regime (initial masses greater than 8 
Mo), we adopt a value of fo, = 0.0415 motivated by the 
work of Brott et al. (2011), who used the step overshoot 
formalism. Both of these values of fj, are measured 
from a distance of 0.008 the local pressure scale height 
into the convection zone from the formal convective- 
radiative boundary. This is the same approach adopted 
in the MIST models (Choi et al. 2016). In order to trans- 
late between the step and exponential-decay versions of 
overshoot mixing, we rely on the work of Claret & Tor- 
res (2017) which shows that the free parameter in the 
step formalism is a factor of ~ 10 larger than fov (their 
Figure 3). For stars with initial masses between 4 Mo 
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and 8 Mo we smoothly ramp between the two values of 
fov. The mass range was chosen to be roughly consistent 
with the ranges considered in the two studies. 

We include no extra mixing due to semiconvection (in 
the sense of Langer et al. 1983), as this process is implic- 
itly accounted for in the convective premixing scheme 
(Paxton et al. 2019). 


4. ADOPTED BINARY-STAR EVOLUTION 
PHYSICS 


POSYDON is predominantly a BPS tool that simulates 
the evolution of an ensemble of binary systems through 
various stages of their life. In the POSYDON framework, 
we base the evolution of binary systems on three ex- 
tended MESA binary grids, as shown in Figure 1. One 
grid consists of initially detached binary systems of two 
H-rich stars starting from ZAMS, where we follow the 
internal evolution of both stars with detailed models 
(Section 5.5). A second grid consists of H-rich stars 
in a semi-detached system with a CO companion (Sec- 
tion 5.6), and a third grid consisting of naked helium 
stars in an initially detached system with a CO com- 
panion (Section 5.7). For the non-CO components in 
these binaries, we follow the same prescriptions for stel- 
lar structure and evolution as described in Section 3. 
However, the internal structures of stars can be affected 
by the presence of a companion, principally through 
tidal interactions and mass transfer. In this section we 
describe how we use the binary module within MESA 
to model each binary's orbit, while self-consistently ac- 
counting for the impact on each star's structure. 


4.1. Tides 


Tidal forces take place in binary systems, as each star 
tends to be deformed by the gravitational pull of its 
companion. Invoked by this gravitational deformation, 
frictional forces inside a star drive a binary toward cir- 
cularization and stellar spin - orbit synchronization. In 
our MESA binary grids, we assume that the initial orbit is 
circularized and the stellar spins are synchronized with 
the orbit. This assumption should be valid especially 
for close orbits of massive stars, where tides are strong 
(Portegies Zwart & Verbunt 1996; Hurley et al. 2002). 
As binaries evolve, their orbital periods change as well as 
the individual stars' rotation periods, potentially driv- 
ing them out of synchronization. Therefore, it is only 
the process of spin-orbit coupling that is relevant for our 
grids of detailed binary-star models. In Section 8.1 we 
discuss our treatment of eccentric, detached binaries. 

We follow the linear approach to tides, which defines 
a timescale for synchronization (Hut 1981). In this ap- 
proach, a torque is applied to non-degenerate stars in a 


binary corresponding to the difference between the or- 
bital and the spin angular velocity AQ, divided by the 
synchronization timescale Tgync: 

AQ 


Tsync 


6Q= 


ÔT, (2) 


where ôr is the timestep and 6 is the change in spin 
angular velocity over a particular timestep. Since every 
layer of a star rotates with its own angular frequency, 
Eq. (2) is separately applied to every layer. The torque 
applied to each layer of the star is added up and an 
opposite torque is applied to the binary’s orbit to en- 
sure angular momentum conservation. Winds and mass 
transfer somewhat complicate the picture, and Paxton 
et al. (2015) gives a detailed description of how these 
effects are accounted for. 

We separately calculate Tsync for both stars (Hut 1981; 
Hurley et al. 2002; Paxton et al. 2015): 


1 kN MF? (R\° 
—=3(F)q I (2) l 9) 


Here M, R and I are the mass, radius and moment 
of inertia of the star for which we calculate the tides, 
q = Meomp/M is the binary mass ratio, and a is the 
orbital separation. k is a dimensionless apsidal motion 
constant characterizing the central condensation of the 
star, and T' is the characteristic timescale for the or- 
bital evolution due to tides. As Eq. (3) shows, Tsync is 
strongly dependent on the ratio of the stellar radius to 
the binary orbital separation. In practice, the quantity 
k/T also varies significantly, depending on whether tidal 
dissipation occurs principally within convective regions 
(due to turbulent friction) or radiative regions (due to 
dynamical tides interacting with stellar oscillations). As 
stars may have both convective and radiative regions 
during their lives, at every timestep taken by MESA we 
separately calculate the dynamical and equilibrium tidal 
timescales, layer by layer, and apply the shorter of the 
two. 

For radiative regions in a star, we calculate k/T' based 
on the dynamical tidal timescale from Zahn (1977, see 
also Hut 1981), where 


k GM R? 
(5) syao A 
rad 


with E> as the second order tidal coefficient and G the 
gravitational constant.^ For the calculation of Es, we 


^ This equation is equivalent to Eq. (42) of Hurley et al. (2002), 
apart from the typo correction of a square root, found in Sepinsky 


et al. (2007). 
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adopt the latest prescriptions from Qin et al. (2018), 
who investigated the dependence of the parameter on 
the convective radius Reony for various metallicities and 
evolutionary stages, finding 


10-942 (Reonv/R)*° for hydrogen-rich stars 
1g- 9 99 B oa E for stripped-helium stars. 


(5) 

For the equilibrium tidal timescale, we calculate the 

synchronization timescale for each convective region in 

the stellar envelope using Eq. (3) and following Hurley 
et al. (2002, cf. Eq. 30): 


(3) zx 2 fconv Mconv.reg. (6) 


E = 


T ^ Qin M ^" 


and use the shortest timescale among them. In Eq. (6), 
Meonv.reg. is the mass of the convective region, feony 
is a non-dimensional numerical factor less than unity 
that takes into account slow convective shells that can- 
not contribute to the tidal viscosity within an orbital 
timescale, and Teony is the convective timescale, which 
we take from Eq. (31) of Hurley et al. (2002, based on 
Rasio et al. 1996), adapted to also accomodate convec- 
tive regions that are below the surface: 


Moony reg. R conv.reg. R conv.reg. 
Teonv =0AS1 aoe) Sense. BOUE 
3L 2 
1/3 
x (Ri conv.reg. = Rb conv.reg.) . (7) 


In the equation above, Riconv.reg.; Rb,conv.reg., and L 
are the radii of the top and bottom boundaries of the 
region, and the stellar luminosity respectively, in solar 
units. We use the surface luminosity in all calculations, 
as it is approximately constant throughout the envelope. 
Typically, the shortest equilibrium tidal timescale cor- 
responds to the outermost convective region. In order 
to avoid tides being dominated by potential artificial 
convective shells that may appear during the numeri- 
cal calculation of a star’s evolution, we only take into 
account regions that consist of at least 10 consecutive 
shells in our models. 

In Figure 2, we show two example evolutionary tracks 
in a Hertzsprung-Russell Diagram of the primary star 
in an interacting binary of Mj initial = 56.46 Mo and 
M2 initial = 28.23 Mo, for two different initial orbital pe- 
riods (3.16 days, top panel; 31.62 days, bottom panel). 
We show which term of the tidal timescale dominates 
the tidal forces: the dynamical tidal timescale with or- 
ange, from Eq. (3) and Eq. (4), assuming the whole 
star is radiative, or the equilibrium tidal timescale with 
blue, from Eq. (3) and Eq. (6), according to the most 


* Core-He depletion 
Equilibrium tides * Core-C depletion 
ZAMS @ RLO 

End of MS 


Dynamical tides 


xx es 
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logio (Tog i /K) 
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Figure 2. Evolution of the donor star in two example 
close binary systems of initially Mijnitiai = 56.46 Mo and 
Mo,nitiai = 28.23 Mo, for two different initial orbital peri- 
ods: 3.16 days (top) and 31.62 days (bottom). The colors 
show which tidal timescale is shortest and dominates in the 
tidal process: equilibrium (blue) or dynamical (orange). Star 
symbols depict the main evolutionary points. 


important convective region of the star. We see that 
in the beginning of the evolution, the dynamical tidal 
timescale dominates, as expected for the MS of mas- 
sive stars that have a radiative envelope. The 3.16 day 
period system (top) initiates early RLO and does not 
form convective layers massive enough for equilibrium 
tides to dominate, until after the end of its MS. For 
the wider binary, even during the MS, the equilibrium 
tidal timescale tends to become comparable to the dy- 
namical tidal timescale, due to convective regions that 
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appear close to the surface of the star. These regions 
include a small part of the total mass of the star (as 
low as 107? Mo), but have a significant radial thickness. 
The equilibrium tidal timescale dominates in all the re- 
maining parts of the evolution, apart from the He core 
burning phase of the stripped primary in the 3 days pe- 
riod system. 


4.2. Mass-transfer 


Over the course of a binary’s evolution, the outermost 
layers of one of the binary’s stellar components may be 
removed, due to the gravitational pull of its compan- 
ion. As a consequence of either the binary’s orbit decay 
or the expansion of a star’s envelope, this transfer of 
mass is dictated by the geometry of the Roche poten- 
tial, namely the gravitational potential constructed in 
the co-rotating reference frame of the binary system. 


4.2.1. Mass Loss Rates from a Star Overflowing its Roche 
Lobe 


To calculate mass-loss rates (due to mass transfer 
only) from main-sequence (MS) stars that overfill their 
Roche lobes, we use the contact scheme within MESA. 
This prescription is a numerical approximation; for stars 
overfilling their Roche lobes at the beginning of each 
timestep, mass is removed such that by the end of the 
timestep, the star remains confined to within its Roche 
lobe. For MS stars, this approximation is consistent 
with more accurate methods, as MS stars are compact, 
with relatively small pressure scale heights. We choose 
this prescription, as it allows us to evolve binary systems 
in which both stars overfill their Roche lobes simultane- 
ously (Marchant et al. 2016). 

As stars evolve off the MS, however, they tend to 
expand, forming less dense envelopes as they become 
giant stars. The large pressure scale heights of giant 
stars cause the contact scheme to become inaccurate, 
and a prescription is required that can more accurately 
treat these stars’ extended envelopes. For stars with a 
central H abundance less than 1079, we switch to the 
Kolb scheme (Kolb & Ritter 1990).° This prescription 
allows the star to expand beyond its Roche lobe, and 
self-consistently calculates the rate that mass can flow 
through the inner Lagrangian point based on the local 
fluid conditions. 


4.2.2. Mass Accretion onto a Non-Degenerate Companion 


5 The current version of MESA does not allow the reversing of the 
donor star in this scheme. Occasionally, an once accreting star 
evolves off the MS, expands, and itself overfills its Roche lobe. 
In these cases mass transfer is not calculated, and the system 
predominantly leads to L2 overflow. 


The evolution of a binary during a mass transfer phase 
depends not only on the mass-losing star but also on the 
mass-gaining star. Based on the nature of the accretor, 
the process of accretion is treated differently. 

For binaries with a non-degenerate accretor (those in 
our grid of two H-rich stars), initially all the mass lost 
by the donor through RLO is accepted by the accre- 
tor. We assume that material being accreted carries the 
specific angular momentum according to de Mink et al. 
(2013, Appendix A.3.3). This prescription allows for 
the distinct treatment of accretion via direct impact of 
the incoming stream on the stellar surface, or, in case 
of accretion onto a more compact star, the formation 
of a Keplerian disk around the accretor. The accreted 
angular momentum spins up the accretor, and mass ac- 
cretion is restricted when the accretor reaches critical 
rotation. Mass falling within a critically rotating ac- 
cretor's gravitational potential will be ejected from the 
binary with the specific angular momentum of the ac- 
cretor (Paxton et al. 2015) in the form of rotationally- 
enhanced stellar winds, following Eq. (1). At the same 
time that accretion is spinning up the outer layers of 
a non-degenerate star, internal mixing processes trans- 
port the surface angular momentum toward deeper lay- 
ers, slowing the star’s rotation rate. 


4.2.3. Accretion onto a Degenerate Companion 


Mass transfer onto a degenerate star proceeds sim- 
ilarly as that onto a non-degenerate star, with a few 
notable exceptions. The primary exception is that mass 
transfer is capped at the Eddington limited rate. For 
sub-Eddington mass-transfer rates onto a CO, mass 
transfer is assumed to be conservative. However, for 
super-Eddington rates, the excess matter is lost from 
the vicinity of the accretor as an isotropic wind (i.e., 
with the specific angular momentum of the accretor). 

We calculate the Eddington-limited rate using stan- 
dard formulae (Frank et al. 2002). We first calculate the 
Eddington luminosity Lgqa for an accretor with mass 
Macc: 

i ene Miast (8) 


K 


where « is the opacity of the incoming material and c 
is the speed of light. For a fully ionized gas, Thomp- 
son scattering dominates the opacity, so x = 0.2(1 + 
X)cm?g~!, where X is the hydrogen abundance of the 
donor. By setting Lgaa equal to the radiation released 
by accreted matter as it falls into a CO's potential well 
(Lace = nM c?), we can recover the Eddington-limited 
accretion rate, 


471 G Mace 
ken 


Mraa = (9) 
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The dimensionless constant 7 sets how efficiently the 
rest mass energy of the incoming matter is converted to 
outgoing radiation, 


2 G Macc 
Te Mat" 


(10) 


For BHs, Race is set by the spin-dependent innermost 
stable circular orbit, while for NSs, we use a constant 
Race of 12.5 km (Most et al. 2018; Miller et al. 2019; Ri- 
ley et al. 2019; Landry et al. 2020; Abbott et al. 2020a; 
Kim et al. 2021; Biswas 2021; Raaijmakers et al. 2021). 
Our grids containing a CO components, are currently fo- 
cused on NS and BH accretors, for which we simulate a 
range of masses. The type of a CO is determined solely 
on its mass, with COs having gravitational mass less 
than 2.5 Mo being classified as NSs, while those with 
mass greater than 2.5 M; as BHs. Within our code, 
the only difference between these types of accretors is 
the corresponding 7; otherwise accretion proceeds iden- 
tically regardless of the type of CO accretor. 

As these COs accrete material, they ought to accrete 
angular momentum. In our current version of the grids, 
we ignore any corresponding increase in the spin rate of 
NSs. For BHs, on the other hand, we self-consistently 
incorporate the increase in spin frequency as well as its 
effect on 7 through the radius of innermost stable circu- 
lar orbit (ISCO; Podsiadlowski et al. 2003), 


Mad 
)—1-—4/1— ( ) , (11) 
where M? 


4cc iS the initial mass of the accreting BH and 
Macc its current one. In this equation it is also implicitly 
assumed that the birth spin of the BH is ~0, as it has 
been suggested by several studies (Fragos & McClintock 
2015; Qin et al. 2018; Fuller & Ma 2019). The cor- 
responding increase in the BH’s non-dimensional spin 
rate, aspin; can be calculated following Thorne (1974) 
and King & Kolb (1999), 


2 1/2 Mi Mi 2 
spin — ace 4 18 2 2 
i " (3) Macc | (ss) 
(12) 


These equations both assume that the spin-up occurs 
due to angular momentum accretion from a disk that is 
truncated at the ISCO. 

We do not explicitly stop our simulations if a NS 
accretes enough mass to cross our 2.5 Mo threshold 
thereby collapsing into a BH, but we do switch the 7 
instantaneously. 


1/2 


4.2.4. Onset of CE 


Stars with radiative envelopes entering RLO respond 
to mass loss by shrinking (Hjellming & Webbink 1987); 
mass transfer then reaches a natural equilibrium set by 
the strength of some driving force (nuclear evolution, 
tides, thermal expansion or some other effect), and how 
quickly the orbital separation, and thus the Roche lobe 
radius, change due to mass transfer through the inner 
Lagrangian point. However, in certain circumstances 
mass transfer increases in a runaway process, either be- 
cause stars expand due to mass loss (e.g., stars with 
deep convective envelopes) or because the binary's orbit 
shrinks faster than a donor star’s radius (e.g., Paczyński 
& Sienkiewicz 1972). These phases of binary evolution 
are notoriously difficult to model as they are intrinsically 
three-dimensional processes, and they span many orders 
of magnitude in spatial and temporal scales (Ivanova 
et al. 2013). We therefore stop our MESA models when 
binaries enter dynamically unstable mass transfer; we 
provide a description of how we address this phase in 
Section 8.2. Here we focus on the conditions we use to 
identify when a binary enters dynamically unstable mass 
transfer. 

First, we assume that a dynamically unstable RLO 
phase is initiated whenever mass-transfer rate exceeds 
0.1 Mc yr-!. It is expected that binaries reaching this 
limit will only further increase their mass-transfer rates, 
as this corresponds to a dynamical limit on the mass-loss 
rate for giant stars (with dynamical timescales of years). 
As a check, we carried out a calibration test where we 
followed the evolution of a test binary to mass transfer 
rates even larger than 0.1 Mo yr^!. In every test we ran, 
we found that the mass-transfer rate increases to arbi- 
trarily high rates, confirming the validity of our limit. 
Assigning a limit to the mass-transfer rate also avoids 
numerical issues caused by the effort of stellar models 
to converge with such extreme mass loss. 

As a second condition, we assume dynamically unsta- 
ble RLO occurs when the stellar radius of the expand- 
ing star extends beyond the gravitational equipotential 
surface, passing through the second Lagrangian point 
(L2). In such cases the lost matter from the L» point 
carries substantial angular momentum, rapidly shrink- 
ing the orbit and leading to a runaway process in which 
the two stars spiral in and trigger a CE (Tylenda et al. 
2011; Nandez et al. 2014). We use the prescription from 
Misra et al. (2020, cf. Eq. 15-19) to define the spherical- 
equivalent radius corresponding to Lə. This condition 
cannot occur for MS donors, since the contact scheme 
for RLO forces a star's radius to be contained within its 
Roche lobe. For cases where two MS stars overfill both 
their Roche lobes, in an over-contact binary, we alterna- 
tively use the prescription from Marchant et al. (2016, 
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cf. Eq. 2) for the Lə radius, which considers that both 
stars can contribute to the overflow of the Lo volume 
together. 

For CO accretors, we set a third condition for unsta- 
ble mass transfer based on the photon trapping radius 
(Begelman 1979; King & Begelman 1999). Inside that 
radius photons are advected inward along with accreted 
matter onto the accretor, while outside that radius, pho- 
tons diffuse away. For stable mass accretion, the photon 
trapping radius occurs close to the accretor; however, as 
the accretion rate increases, the photon trapping radius 
expands. Once the photon trapping radius reaches the 
Roche-lobe radius of the accretor it is assumed to lead 
to a CE phase. Since the radius of the photon trap- 
ping envelope Rtrap depends on the Eddington limit of 
the accretor Mgaa and the mass-transfer rate from the 
donor Maäonor, we limit the latter assuming an instability 
condition when (Begelman 1979): 


Maonor > Mpg =e, (13) 


As a final condition, occasionally two stars will both 
overfill their Roche lobe while one of those stars has 
evolved off the MS. Since the contact scheme in MESA 
can only evolve binaries in which both stars are on the 
MS, we assume these binaries automatically enter a CE. 

As a test we compared the first three dynamical in- 
stability conditions separately to investigate their effect 
and found that the limiting mass accretion rates are all 
similar: as soon as a binary reaches any one of them, the 
other two are close to their limits as well. Therefore, for 
a particular binary in our MESA simulation, the binary 
is considered to enter a CE if any one of them occurs. 


5. GRIDS OF DETAILED SINGLE- AND 
BINARY-STAR EVOLUTION MODELS 


While in Section 3 and Section 4 we describe the 
physics we adopt in our simulations, here we provide 
numerical details about how we produce each of our five 
MESA grids. This includes our procedure for producing 
initial stellar models for each of our grids (Section 5.1), 
our termination conditions common across all of our 
grids (Section 5.2), and a description for each of our 
five grids of binary simulations (Sections 5.3-5.7). We 
summarize the basic properties of each grid in Table 2. 


5.1. Zero Age Main Sequence Models 


We create our own library of ZAMS models for both 
H-rich and He-rich stars. For the creation of the H-rich 
ZAMS models we use the MESA revision 11701 template 
create_zams. The process begins with creating a fully- 
convective star with no nuclear fusion taking place and 


adopting our protosolar abundances (Section 3.1). This 
model is then evolved with our adopted nuclear reaction 
network until the H-burning luminosity exceeds 99% of 
the total luminosity. 

The He-rich ZAMS (ZAHeMS) models are created in 
three steps. First, we create a pre-MS He star with 100% 
tHe in the same way that we create a H-rich pre-MS star. 
In the second step we adjust the initial metallicity. In 
the third step we evolve the model until the He-burning 
luminosity exceeds 99% of the total luminosity. 

The two sets of ZAMS and ZAHeMS models are used 
as a starting point for the five grids of single- and binary- 
star models. 


5.2. Termination conditions 


We set conditions for the termination of our evolu- 
tionary models based on both single-star properties and 
binary-star properties. If any one of these conditions 
are met by an individual simulation, it is terminated at 
that timestep. Our termination conditions are: 


e A star’s age exceeds the age of the Universe 
(13.8 Gyr), a condition that is typically only met 
for the lowest-mass stars we simulate (Minit SS 
0.8 Mọ). In our single-star grids, for numerical 
purposes, we allow stars to evolve beyond this con- 
dition, then truncate their evolution afterwards at 
the end of the MS, which may extend beyond the 


age of the Universe. 


e A star becomes a WD, a condition we quantify by 
checking if the central degeneracy parameter Te 
(Coulomb coupling parameter) exceeds 10 (Choi 
et al. 2016). 


e A star reaches the end of core C-burning, a con- 
dition triggered when the fractional abundances 
of both C and He decrease below 107? and 10-6, 
respectively, at the star's center. 


e A binary enters a CE phase, as described in Sec- 
tion 8.2. 


e A star reaches the thermally-pulsating asymptotic 
giant branch (TP-AGB) phase and then reaches a 
point of failure (numerical non-convergence) dur- 
ing thermal pulsations. Because these are not un- 
common and we consider the nascent WD to be 
well-formed be within the AGB star, we consider 
this evolution to be successful. 


Our simulations may occasionally end prematurely be- 
fore any of the aforementioned conditions are reached. 
'This may happen because the minimum timestep limit 
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Table 2. Summary of the five detailed single- and binary-star model grids. 


Initial state Parameters’ range and resolution 


Star 1 Star2 Mi[Mgo] Alogyyg Mi M2[Mo] Alogi 9 M2 q Aq Porb [day] A logio Porb NA Failures 2 
ZAMS - 0.5—300 0.014 - - - - - - 200 1.596 
ZAHeMS^ - 0.5-80 0.055 - - - - - - 40 096 
ZAMS ZAMS 76.2-120 0.025 - - 0.05-0.95 0.05 0.72-6105 0.07 56000 58240 20-9 1.596 
Evolved, H-rich d CO 0.5-120 0.06 1-35.88 0.074 - - 1.26-3162 0.13 25200 8-8 0.9% 
ZAHeMS CO 0.5-80 0.055 1-35.88 0.074 - -  0.02-1117.2 0.09 39480 47 4.8% 


a Total number of models in this grid. 


b Percentage of models that stopped due to numerical-convergence errors before reaching one of our stopping conditions. These rates describe the 
finalized grids, after a series of re-runs have occurred; see Section 6.1 for details. 


© Zero-age He Main Sequence stars. 


d Although this grid is initialized with H-rich stars at ZAMS, we ignore the portion of each simulated binary's evolution prior to the onset of RLO. 
'The initial state of Star 1 in this grid are therefore somewhat evolved. 


within MESA (10^? s) is reached or any individual simu- 
lation reaches our maximum runtime on our computing 
cluster (set to 48 hours). We provide details describing 
how we approach such runs in Section 6.1, but these 
failures are rare, occurring a few percent or less in each 
grid. 

All single- and binary-star models in POSYDON have 
a final, internal structutre profile written to correspond 
precisely to the last evolutionary phase milestone. These 
profiles are discussed further in Section 6.2, Section 8.2, 
and Section 8.3. 


5.3. H-rich, single-star grid 


Our first grid of single-star evolutionary models con- 
tains a series of non-rotating H-stars with our adopted 
protosolar composition of Y = 0.2703 and Z = 0.0142. 
The grid consists of 200 masses, ranging from Minit = 
0.5 Mo to Minit = 300 Mo with a logarithmic spacing of 
Alogig(Minit/Mo) = 0.014dex. For each star, models 
were initialized using the procedure described in Sec- 
tion 5.1, and evolved until one of the termination con- 
ditions provided in Section 5.2 occurs. 

To test their validity, we compare POSYDON evolution- 
ary tracks to the widely used stellar evolution tracks 
from the Geneva (upper panel; Ekstróm et al. 2012), 
MIST library (center panel; Choi et al. 2016), and BSE as 
implemented in COSMIC (lower panel; Pols et al. 1998; 
Hurley et al. 2000; Breivik et al. 2020) groups in Fig- 
ure 3. In all cases we show non-rotating models with 
initial masses between 1 Mọ and 300 Mo. The pre-MS 
evolution is omitted from the MIST evolutionary tracks 
and the TP-AGB and post-AGB phases are omitted for 
clarity. All sets of tracks show a similar location for 
the ZAMS; the subsequent evolution along the MS and 


through He-burning phases differs due to the way each 
set of models treats mixing across convective bound- 
aries. The clearest differences between the POSYDON and 
other models are in the location of the hook feature near 
the MS turnoff for higher masses, a result of the differ- 
ent adopted core overshoot treatments, and the posi- 
tions of later phases, a result of the different wind mass- 
loss treatments among the different groups. The COSMIC 
tracks extend to larger radii and cooler effective temper- 
atures, which may place them in the regime of LBVs; 
however, none of the other sets of evolutionary tracks 
enter this regime. For a more in depth comparison see, 
e.g., Agrawal et al. (2020, 2022). 

Figure 4 compares the final C/O core mass between 
POSYDON, MIST, and SSE as implemented by COSMIC. 
SSE models are calculated until central C-burning, while 
MIST and POSYDON are calculated through central C- 
exhaustion for those stars with sufficient mass to ig- 
nite carbon or to the WD cooling sequence for lower 
masses. Differences between the core masses of MIST and 
POSYDON are generally due to the different overshooting 
parameter (we adopt foy = 0.0415 for stars with masses 
above 8 Mo, compared with fo, = 0.016 adopted by 
MIST). The COSMIC/SSE models exhibit different behav- 
ior at larger masses as these prescriptions are based on 
stellar models that only go up to 50 Mo; larger masses 
than this are an extrapolation. 

As a final comparison, Figure 5 shows a Hertzsprung- 
Russell diagram of a subsample of POSYDON single stel- 
lar model where we indicate the different evolutionary 
POSYDON stellar states (Figure 17) across a range of stel- 
lar masses. 


5.4. He-rich, single-star grid 
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Figure 3. Comparing POSYDON H-rich ZAMS evolutionary 
tracks (blue in all panels) with non-rotating Geneva (upper; 
Ekstróm et al. 2012), MIST (middle; Choi et al. 2016), and SSE 
with default wind mass loss prescription from COSMIC version 
3.4.0 (lower; Pols et al. 1998; Hurley et al. 2000; Breivik et al. 
2020) tracks in the Hertzsprung-Russell diagram. The mass 
range shown is 1-300 Mo in all cases. The masses shown are 
the same as the Geneva grid of models between 1 and 120 Mo 
with the addition of 175 and 300 Mo for POSYDON, MIST, and 
COSMIC/SSE. The dotted, gray lines indicate constant radius 
at powers of 10 in Ro. 


Our second grid of single-star evolutionary models 
consists of non-rotating He-rich stars with Y; = 
1 — ŽZinit and our adopted protosolar Zi, = 0.0142. 
This grid consists of 40 masses ranging from Minit = 
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Figure 4. Comparison of the final C/O core mass in 
SSE as implemented by COSMIC (magenta), MIST (green), 
and POSYDON (blue). Differences between MIST and POSYDON 
are due to the larger core overshoot parameter adopted by 
POSYDON. Disagreement with the SSE models is expected as 
these models are based on simulations that were only com- 
puted for initial masses up to 50 Mo; more massive stars are 
an extrapolation. 


0.5 Mo to Minit = 80 Mo with a logarithmic spacing of 
A logi9(Minit/Mo) = 0.055 dex. For these masses, stel- 
lar evolution models were computed starting from ZA- 
HeMS models (Section 5.1) and evolved until one of the 
termination conditions described in Section 5.2 occurs. 
For all but the lowest-mass cases, the core C depletion 
condition is the relevant one; models with initial masses 
below 1.1 Mc do not ignite C-burning in the core, and 
therefore terminate as He-core WDs. 

As a test of our POSYDON He-star models, we com- 
pare their lifetimes to those of the He-star models of 
Woosley (2019) in Figure 6. Only the overlapping range 
of initial masses is shown here; the Woosley (2019) grid 
includes models with masses from 1.8-120 Mọ. These 
models match to within ~ 0.1 dex in log lifetime across 
the entire range of initial masses. 

As a second test, we compare the final masses between 
the same two model grids in Figure 7. Although the 
lifetimes are similar, the final masses show a significant 
difference, particularly at higher initial He-star masses. 
Woosley (2019) notes that the change in slope of the 
initial-final mass relation around Mgya4| of 11 Mọ is due 
to the mass-loss prescription adopted for exposed CO 
cores; at larger initial He-star masses, the entire He-star 
mass is burned to heavier elements. For all the single- 
and binary-star model grids in POSYDON, we adopt the 
mass-loss prescription from Nugis & Lamers (2000) for 
He-rich stars. The latter predicts on average stronger 
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Figure 5. Hertzsprung-Russell diagram of a subsample of 
POSYDON single-stellar models where the different POSYDON 
stellar states are indicated according to the legend. 
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Figure 6. Lifetimes of the POSYDON and Woosley (2019) 
single He-star evolution models match to within ~0.1 dex. 
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Figure 7. Final masses of the POSYDON and Woosley (2019) 
single He-star evolution models. Differences at Miniti Z 


N 


20 Mo are due to the stronger wind mass loss prescription 
adopted by POSYDON. 


Figure 8. Evolution of radius for He stars with masses 
between 2 and 10 Mo. Less massive He-stars evolve slower, 
but expand farther when they become giant stars. This trend 
implies that more-massive He-stars in binary systems will 
undergo RLO for a relatively narrow range of orbital periods, 
a behavior exhibited by our binary star simulations and seen 
in Figure 14. 


wind mass loss than the prescription from Yoon (2017) 
adopted by Woosley (2019), leading to the substantially 
different final masses between the two prescriptions at 
Minitial Z 20 Mo. 

Finally, we show the radius evolution of the He star 
tracks for the mass interval of 2-10 Mo in Figure 8. He- 
rich stars exhibit a peculiar feature where less-massive 
stars expand farther on the giant branch than their more 
massive counterparts, in agreement with results by Ha- 
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bets (1986). When in a binary system, this implies that 
there is a relatively narrow range of orbital periods in 
which massive He stars will undergo RLO. Less massive 
He stars expand to hundreds of Ro, leading to a wide 
range of orbital periods in which these stars can interact 
with a putative companion. This behavior is realized in 
our binary star grids, and its effects are seen explicitly 
in Figure 14. 


5.5. Binaries consisting of two hydrogen-rich 
main-sequence stars 


For modeling the evolution of two ZAMS stars in 
a binary system, we run a grid of 58,240 separate 
binary evolution models, varying the initial mass of 
the primary star Mı, the initial binary mass ratio 
q = M2/Mı (where Mə is the mass of the compan- 
ion star), and orbital period Psr». We consider 52 
values of initial primary masses, ranging from Mı = 
6.22 Mo to Mı = 120 Mọ with a logarithmic spac- 
ing of Alog,)(Mi/Mo) = 0.025 dex, and 20 values of 
initial binary mass ratios, ranging from q = 0.05 to 
q = 1 with a spacing of Aq = 0.05. Finally, we 
cover 56 values of initial orbital period, ranging from 
Pay = 0.7 days to P4, = 6105 days with a logarithmic 
spacing of Alog,9(Porn/days) = 0.07 dex, in order to 
explore all binary configurations ranging from close sys- 
tems in initial RLO to wide systems that never exchange 
any mass. 

We simulate binaries by first separately initializing 
two H-rich, single stars at ZAMS following the proce- 
dure defined in Section 5.1. We then place those stars 
in a binary with a second relaxation step, where we force 
their their rotation periods to be synchronized with the 
orbital period, implicitly assuming that the synchroniza- 
tion has happened during the pre-MS phase. The latter 
might not be true for wide binaries, but our assumption 
induces negligible rotation to the stellar components of 
those systems and does not affect their further evolution. 
As long as both stars in the binary are under-filling their 
Roche lobes after this relaxation step, we start to evolve 
the binary. Evolution continues until one of the termi- 
nation conditions described in Section 5.2 occurs. 

In Figure 9 we provide two two-dimensional slices of 
this grid, where we show our simulation outcomes as a 
function of Mı and Pa for fixed q values. In the left 
panel we show one example of a mass ratio q — 0.3, and 
on the right a more-equal mass q — 0.7 slice. Each point 
in the panels represents a separate simulation from our 
grid. Diamond markers represent models that termi- 
nated in a CE, while square markers represent models 
that terminated when one of the stars completed its evo- 
lution (e.g., reached core C exhaustion). These are sys- 


tems that experienced either only stable mass-transfer 
episodes or no mass transfer at all, so their evolution can 
be continuously modeled. At the bottom of each panel 
we see systems that are born filling their Roche lobes 
(black dots). These systems are assumed to merge, and 
therefore never produce a viable binary. Finally, a small 
fraction of systems never complete their evolution, pro- 
ducing binary stellar models that at some point fail to 
converge (red diamonds). 

Separately, the color of each marker indicates that par- 
ticular binary's mass transfer history. Systems with suf- 
ficiently close initial Porp tend to lead to contact phases 
(orange) where both stars fill their Roche lobes simul- 
taneously. Most, but not all, of these system end up 
entering a CE phase. Sufficiently widely separated (or 
very massive) systems never fill their Roche lobes, and 
therefore never interact (gray markers). For interme- 
diate orbital periods, the colors differentiate the evolu- 
tionary state of the donor when the latest mass transfer 
phase was initiated, ranging from MS (blue) to post-MS 
(tan), to stripped He-MS (brown). Stable mass transfer 
causes the donor star to be almost completely stripped 
of its H-rich envelope. In the latter case (brown) the 
low-mass stripped donors initiate a second mass trans- 
fer phase (Case BB mass-transfer) when they re-expand 
(Delgado & Thomas 1981; Laplace et al. 2020). 

Comparison between the two panels shows that the 
mass ratio leads to a stark difference in the mass-transfer 
outcomes. Whereas nearly all systems with an initial 
q = 0.7 result in stable mass transfer, the opposite is 
true for our q = 0.3 systems. At the same time, some 
features between the two mass ratios are similar: (i) 
The boundary between interacting and non-interacting 
systems seems to be insensitive to q (and therefore the 
secondary's mass). At the largest orbital periods, stars 
do not expand far enough to overfill their Roche lobes. 
At the largest masses, stars have extremely strong winds 
that widen their orbits, simultaneously stripping the pri- 
mary of its H-rich envelope, and these stars never ex- 
pand enough to fill their Roche lobes. (ii) Systems with 
initial Porp S; 5 days tend to result in dynamically un- 
stable mass transfer. (iii) There is a large region of bi- 
naries with initial primary mass 740-50 Mo that stably 
overfill their Roche lobes as post-MS stars. These stars 
achieve their mass transfer stability mainly due to their 
strong stellar winds, which increases the mass ratio and 
the orbit of the system until the moment of overflow. 

We model, and keep track of, the properties of both 
stars in the binary system throughout their evolution, as 
well as their detailed internal structure at the end of the 
models. In Figure 10 we show, for the same two mass- 
ratio slices as in Figure 9, the final rotational rate of the 
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View of two grid slices, for two different values of initial binary mass ratio (q = 0.3 on the left, q = 0.7 on the 


right), from our grid of binary-star models consisting of two H-rich stars, initially at ZAMS. The different symbols summarize 
the evolution of each of the models. We distinguish between models that experienced stable or no mass transfer (squares), 
reaching the end of the life of one of the stars, and the ones that stopped during mass transfer due to one of our conditions for 
dynamical instability (diamonds). Different colors distinguish the evolutionary phase of the donor star during the latest episode 
of mass transfer (or no RLO at all for grey). Small black dots at low initial periods depict systems that were in initial RLO at 
birth and red diamonds represent the models that stopped prematurely for numerical reasons. 


secondary (the initially less massive) stars for systems 
that avoid dynamically unstable mass transfer. Each 
marker’s color is set by how close each star’s rotation 
rate is to its critical rate. Highly rotating secondary 
stars have all experienced substantial mass and angular- 
momentum accretion during their evolution. Many of 
them have reached critical rotation, (ws/ws,crit)2 = 1, 
early during mass transfer, at which point further mass 
accretion becomes non-conservative (c.f. Section 4.2.2). 
The right-hand panel shows that the companion’s rota- 
tion rate is closely linked with Mj, as companion stars 
with lower mass primary stars also have lower masses 
and therefore do not lose as much angular momentum 
through their own stellar winds. This behavior is inde- 
pendent of the assumed initial rotation of the stars. 
We find a small subset of initially very close systems 
in the bottom right corner (logıo(Mı/Mọo) > 1.75 and 
logio(Porb/days) ~ 0.5) that retain a significant rota- 
tional rate even though they avoid mass transfer. In 
binaries with such tight orbits, tidal forces between the 


stars are sufficiently strong to keep them fast rotating, 
despite their strong winds. 


5.6. Binaries consisting of a compact object and a 
hydrogen-rich star, at the onset of Roche-lobe 
overflow 


Our second grid of binary star simulations consists of a 
H-rich star in a binary with a CO at the onset of RLO. 
This grid consists of 25,200 binary evolution models, 
where we vary the initial mass of the primary star Mı, 
the initial mass of the CO Moo, and the orbital period 
Pos. We consider 40 values of initial primary masses, 
ranging from Mı = 0.5 Mo to Mı = 120 Mọ with a 
logarithmic spacing of A log; )(M1/Mo) = 0.06 dex, and 
21 values of initial CO masses, ranging from Mco = 
1 Mo to Mco = 35.88 Me with a logarithmic spacing 
of Alog;4(Mco/Mg) = 0.074dex. Finally, we cover 
30 values of initial orbital period, ranging from P5, = 
1.26 days to Pop = 3162 days with a logarithmic spacing 
of Alog,,(Porp/days) = 0.13dex. Our choice of CO 
mass range covers massive WD, NS, and BH accretors. 
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rate, (ws/ws,crit)2. In most cases where mass transfer occurred, the secondary star accreted mass and spun up, remaining highly 


we depict the final ratio of the angular velocity of the secondary star (the initially less massive) divided by its critical rotation 
spinning until the end of the life of the initially more massive star. 


Figure 10. For the same grid slices shown in Figure 9, and only for systems where one of the two stars reached the end of its life, 


RLO are not considered further; these detached bina- 
ries are modeled separately as described in Section 8.1. 
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Figure 11. View of two slices, for two different values of initial CO masses (Mco = 1.43 Mo on the left, Mco = 14.66 Mo on 
the right), from our grid of binary-star models consisting of H-rich star and a CO at the onset of RLO. The different symbols 
summarize the evolution of each of the models, as in Figure 9. Binaries that never initiated mass transfer are not shown here. 


stars have winds too strong to expand into giant phases. 
In this grid, Figure 11 shows an additional region of 
white space at low mass (M € 1Mo) that occurs be- 
cause these stars remain on the MS for the entirety of 
the simulation, never expanding to fill their Roche lobes 
within the age of the Universe. 

Examining the stability of the mass-transfer phase, 
Figure 11 shows that nearly every donor star accret- 
ing onto a 14.66 Mọ BH does so stably, whereas only 
the lower mass accretors (M <4.5Mo) do so for NS 
accretors. This difference is because the stability of a 
mass transfer in a binary primarily depends on the mass 
ratio, with a higher accretor mass allowing for higher 
donor masses. Our findings, at least for the case of 
NSs, are consistent with recent results from Misra et al. 
(2020), who use the same criteria to define the onset of 
La overflow leading to dynamical instability as done in 
our work. 

Figure 12 shows the relative changes in the accretor 
masses in the same two slices in Mco as Figure 14. High 
amounts of accretion mainly depends on two factors: a 
sufficiently high-mass accretion rate and a long-lasting 
RLO phase. In both panels, this happens for binaries 
with short periods ~ 1 day, and pre-RLO mass ratios in 
the range q ~ 1-2 (defining q = Mco;/Mi;). Despite 


our assumption of Eddington-limited accretion, for these 
binaries, stable accretion occurs for over a long time, 
and in both cases the binaries transition to low-mass 
X-ray binaries. These findings are in agreement with 
earlier works by Podsiadlowski et al. (2003); Fragos & 
McClintock (2015); Misra et al. (2020). 

'The high mass-transfer rates achieved by most initial 
binary configurations are explicitly shown in Figure 13, 
where each marker's color corresponds to the peak mass- 
transfer rate for each binary. These rates refer to the 
mass being lost by the donor star due to RLO; accretion 
onto the accretor is still Eddington-limited. In both pan- 
els, super-Eddington mass-transfer rates occur in most 
binaries with higher peak mass-transfer rates encoun- 
ters in binaries with higher periods and larger donor star 
masses. However, since the larger orbital separation of 
these binaries implies the donors in these systems would 
be more evolved at RLO onset, compared with initially 
shorter-period binaries, these mass transfer phases tend 
to be short-lived. Therefore, binaries with short orbital 
periods (but not so short that they overfill their Roche 
lobes initially) will lead to the most accretion onto a 
CO. 
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Figure 12. Relative increase in the mass of the CO (Mco, — Mco,i)/Mco, due to accretion for systems where the non- 
degenerate star reached the end of its life. The grid slices are the same as shown in Figure 11. Although accretion is Eddington- 
limited, COs in binaries with pre-RLO mass ratios in the range q ~ 1-2 (defining q = Mco,;/Mi,;) and short initial periods, 
which will experience long-duration mass-transfer phases, manage to accrete a significant amount of mass. 


5.7. Binaries consisting of a compact object and a 
He-rich star 


Our final grid of detailed binary-star simulations con- 
sists of 39,480 models of He-rich stars with CO compan- 
ions, where we vary the initial mass of the primary star 
Mı, the initial mass of the CO Mco, and the orbital 
period Porp. We consider 40 values of initial primary 
masses, ranging from M; = 0.5 Mọ to Mı = 80 Mo with 
a logarithmic spacing of Alog;,(Mi1/Mg) = 0.055 dex, 
and 21 values of initial CO masses, ranging from Mco — 
1 Mọ to Mco = 35.88 Mọ with a logarithmic spac- 
ing of Alogj)(Mco/Mo) = 0.074dex. Finally, we 
cover 47 values of initial orbital period, ranging from 
Pay = 0.02 days to P5, = 1117.2 days with a logarith- 
mic spacing of A log,,(Porp/days) = 0.09 dex. Our pro- 
cedure for generating these binaries closely follows the 
process described in Section 5.5 for the grid of binary- 
star modes composed of two H-rich stars. Here, we re- 
place the initial primary-star models with He-rich stars 
at ZAHeMS, while the companion COs are modeled as 
point masses. 


Figure 14 shows an example of two slices of this grid, 
one corresponding to a NS companion (with Mco œ 
1.43 Mo ) and one corresponding to a BH (with a Mco ~ 
14.66 Mo). Marker shapes and color scheme follow the 
same convention as in Figure 11, but since these simu- 
lations are initialized with He-stars, the symbol key is 
simplified in Figure 14. 

When comparing the two panels, the most appar- 
ent difference occurs at large M; and short orbital pe- 
riod: whereas accreting NSs enter unstable mass trans- 
fer (these systems typically end up merging in a CE, cf. 
Section 10), the corresponding accreting BHs typically 
either overfill their Roche lobes at ZAMS or avoid mass 
transfer altogether. In contrast, we find that indepen- 
dently of the CO mass, systems with low He-star masses 
(Mı € 3 Mo) mass transfer up to wide orbital periods 
(Py < 10? days). This occurs because low-mass He- 
stars expand their He-rich envelope much farther during 
their later He-shell and C-burning phases (Figure 8). 

Both slices of the grid present two islands of failed 
simulations, one with M, c 1.8 Mọ and P4 of the 
order of days and another island with Mı < 1 Mo and 
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Figure 13. Same as Figure 12, but now the color of the symbols depict the maximum mass-transfer rate that occurred in the 
evolution of each binary. A significant part of the parameter space leads to highly super-Eddington mass-transfer rates, albeit in 
short-lived phases for most cases, and thus to the potential formation of ultra-luminous X-ray sources. This peak mass-transfer 
rate refers to the rate the donor star is losing mass through RLO; accretion onto the accretor is limited to the Eddington rate. 


Por» of the order of hours. MESA has difficulty modeling 
the envelope's structure as it expands to large radii in 
the first island, whereas the second, short- P,,4, island is 
due to MESA having difficulty following a star's evolution 
into a He WD after it has been spun up due to tides 
and mass-transfer. Combined, failed runs account for 
c 596 of the models in this grid. In practice we find 
these failed runs do not bias our population synthesis 
results of merging NSs and BHs as these portions of the 
parameter space predominantly lead to the formation of 
WDs. 

In Figure 15 we show the same two grid slices, but 
now the marker color corresponds to the specific angular 
momentum of the He-star j1, at the end of the simula- 
tion. MESA allows us to track this quantity, as it self- 
consistently models the interplay between tides (which 
spin up the star), stellar winds (which spin down the star 
and widen the binary), mass transfer (which alters the 
orbital period), and internal angular momentum trans- 
port. Comparing Figure 14 and Figure 15, we find that 


the He-stars with the highest specific angular momenta 
are those with either short P r» or stable mass transfer. 
'The binary-star grid, composed of a He-rich star and 
a CO companion, presented in this section closely agree 
with those of Qin et al. (2018) and Bavera et al. (2020, 
2021). In contrast to these previous works, the present 
grid further expands the parameter space coverage to 
lower He-star masses and to larger orbital periods. 


6. GRID POST-PROCESSING 


Each single- or binary-star evolution simulation pro- 
duces a series of data files which must be parsed, an- 
alyzed, and collated before we can use them within 
POSYDON. Our process includes: (1) re-running any failed 
simulations; (2) adding post-processed quantities to our 
data grids; (3) a post-processing procedure used exclu- 
sively on our single, H-rich and He-rich star grids, which 
allows for an efficient interpolation among tracks of dif- 
ferent masses; (4) the downsampling of our grids to re- 
duce data size; (5) classifying each model within our 
grids based on the different resulting stellar and binary 
types, and (6) fitting classifiers and interpolators over 
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Figure 14. View of two grid slices for two different values of initial CO masses (Mco = 1.43 Mo on the left, Mco = 14.66 Mc 
on the right), from our grid of binary-star models consisting of He-rich stars and CO companions. 'The different symbols 


summarize the evolution of each of the models, as in Figure 9. 


the stellar and binary parameters in each grid. We de- 
scribe the first 4 steps next, while the steps of classifi- 
cation and interpolation are discussed in Section 7. 


6.1. Re-running Failed Models 


After having computed our grids of single- and binary- 
star models, we first identify those runs that did not 
reach our desired end point (cf. Section 5.2). This can 
happen for a variety of reasons, many of which we have 
not yet been able to eliminate. For example, one source 
of problematic runs appears to deal with stellar oscil- 
lations; in certain cases, MESA tries to resolve short- 
timescale evolution driven by the &-mechanism, which 
dramatically shortens the size of successive MESA steps. 
We address this problem by re-running our failed binary 
simulations with a maximum radiative opacity (Kmax) 
set to 0.5 cm? g~t. This approximation reduces the fail- 
ure rate of each binary grid from ~ 10.996, ~ 8.0%, and 
~ 11.8%, for the binary-star grids composed of two H- 
rich stars, a H-rich star with a CO companion at the 
onset of RLO, and a He-rich star with a CO compan- 
ion, respectively, to ~0.9%, ~1.5%, and ~4.8%. The 
differences in the resulting evolutionary tracks with and 
without the opacity limit are generally small when com- 
pared to differences in tracks of adjacent points in our 


initial parameter space and compared to our interpola- 
tion accuracy (Section 7.5). 

Figure 16 shows a typical example of a binary-star 
model, initially composed of two H-rich ZAMS stars 
with masses 10.50 Mọ and 5.25 Mọ and an orbital pe- 
riod of 43.94 days. This binary initially failed to reach 
the end of the simulation (dashed, black line; MESA ex- 
ceeded its minimum timestep limit), but did so success- 
fully when re-run with an upper limit to the radiative 
opacity (orange line). The top panel shows that the stel- 
lar radius evolves similarly between the two simulations 
as the donor star loses mass. For the radius and ef- 
fective temperature (bottom panel), the two properties 
most affected by an opacity limit, differences between 
the two tracks are typically less than 0.1 dex. 

For further comparison we show an adjacent binary- 
star model in the same grid with stars of the same mass 
ratio and orbital period but slightly less massive primary 
star of Mı = 9.9 Mo, that successfully reached the end 
of the simulation without the need to limit the opacity. 
Although the neighboring simulation is better able to 
match our failed run when the luminosity dips to low 
values, comparison between all three tracks in Figure 16 
suggests that any inaccuracies accrued by our opacity 
limit are of a similar magnitude to any differences be- 
tween adjacent simulations in our model grids. 
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Figure 15. The final specific angular momentum jı = Jı/Mı where Jı is the He-star AM and Mi its mass, at carbon depletion 
for our grid of He-stars with CO companions. We only show jı for systems where the non-degenerate star reached the end of 


its life. The grid slices are the same as shown in Figure 14. 


Table 3. A list of the different stellar types we adopt in POSYDON. 


Compact Non-degenerate star states 

WD H-rich, Core H burning H-rich, Core He. burning H-rich, Shell. H. burning 

NS H-rich, Central, He depleted H-rich Central, C depletion H-rich, non burning 

BH stripped. He Core H burning Stripped. He. Core He. burning stripped. He Shell H burning 


Stripped He Central He, depleted 


stripped He Central. C depletion 


stripped. He. non burning 


6.2. Post-Processed Quantities 


Our single- and binary-star MESA simulations result in 
two types of files: history files which contain the time 
evolution of the binary's and its component stars! prop- 
erties and profile files which define each non-degenerate 
star's structure. In this first version of POSYDON, we save 
the profiles of stars only at the end of the simulations. 
Combined with the MESA terminal output, we have all 
the information necessary to analyze each simulation. 

As a first step, we analyze the final binary proper- 
ties and terminal output to broadly determine how and 
why each binary simulation ended. We first identify the 
small subset of binaries in which, despite the process 


described in Section 6.1, the MESA simulation failed to 
converge; these are ignored throughout the remainder 
of this work. For the successful binary-star simulations, 
we have four separate conditions (e.g., Figure 15): (1) 
binaries that were already in RLO when initialized at 
ZA(He)MS; (2) binaries in which one star reached the 
end of its lifetime (i.e., one of the first three termination 
conditions described in Section 5.2) and went through 
stable mass transfer; (3) binaries which avoided mass 
transfer and the two component stars essentially evolved 
in isolation; and (4) binaries that entered unstable mass 
transfer (described in Section 8.2). 
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Figure 16. Typical example of the evolution of a binary- 
star model that failed to reach the end of the simulation due 
to over-resolved stellar oscillations that eventually lead MESA 
to convergence problems. We show the evolution of the pri- 
mary’s radius, as a function of its mass (top panel) and its 
track in the Hertzsprung-Russell diagram (bottom panel). 
The binary initially consists of a 10.50 Mọ and a 5.25 Mọ 
H-rich ZAMS stars at an orbital period of 43.94 days. Com- 
parison between our original, failed simulation (dashed, black 
line) and a successful simulation where we artificially limit 
the radiative opacity to 0.5 cm?/g (orange line) shows that 
our approximation is typically accurate to within 0.1 dex 
for all stellar parameters. We additionally show the evolu- 
tion of an adjacent simulation in our binary grid (blue line; 
same mass ratio and orbital period but a primary mass of 
9.91 Mo), which shows that any inaccuracy induced by our 
opacity approximation is a similar magnitude to the differ- 
ences between neighboring simulations. 


As a second step, we separately analyze each star and 
assign it a stellar type at the end of each simulation. 
For COs, this is a straightforward task, as the type of 
CO is dependent only on its mass. Non-degenerate stars 
present a more difficult object to typecast. Most pBPS 


codes rely upon the k-type stellar classification intro- 
duced in Hurley et al. (2000, based on Tout et al. 1997) 
In POSYDON, we use instead a two-term classification sys- 
tem, dependent upon both what part of the star (if any) 
is undergoing nuclear burning and the envelope's com- 
position. Table 3 provides a list of the possible stellar 
type combinations, while in Figure 17 we show the al- 
gorithm for determining them. 

As a third step, we assign a designation to the re- 
sulting binary configuration in each of our simulations. 
'These include detached for binaries in which both stars 
are confined within their respective Roche lobes, RLO1 
or RLO2 for binaries in which the primary or secondary 
star, respectively, is overfilling its Roche lobe; contact 
for binaries in which both stars are overfilling their re- 
spective Roche lobes; not. converged for systems where 
the binary-star simulations ran into numerical conver- 
gence problems, and initial.MT for binaries initialized 
in RLO. These designations of stars and binary states 
are used throughout POSYDON and are updated by each 
evolutionary step. Therefore evolutionary phases mod- 
eled with on-the-fly calculations (see Section 8) also af- 
fect the star and binary states, which result in two ad- 
ditional possible designations: merger for those binaries 
that have merged, and disrupted for those binaries that 
have become unbound due to some process. 

As a fourth step, we analyze the mass-transfer his- 
tory of the binaries we simulate, identifying the donor 
star's state when RLO initiated, and whether or not 
that mass transfer phase was stable or unstable (e.g., 
caseA from starl). Note that we use the canonical def- 
initions for Case A, Case B, and Case C mass transfer 
(for a review, see e.g., Iben 1991). Specifically these 
labels (which can be identified by the differing sym- 
bol types and colors shown in, e.g., Figure 11) identify 
whether the donor star was on the MS, on the post-MS, 
or a stripped He star. In cases where systems evolve 
through multiple phases of mass transfer, all phases 
are included in the label (e.g., caseA/B from stariifa 
Case A mass-transfer phase is followed by Case B one). 

As a final step, we calculate a number of post- 
processed quantities, ranging from parameters related to 
specific core-collapse mechanisms (Section 8.3) to differ- 
ent CE prescriptions (Section 8.2). These are typically 
parameters that require integrals over all or part of a 
star's structure, which for efficiency we pre-compute. All 
post-processed quantities are summarized in Table 8. 


6.3. Resampling of single-star grids using Equivalent 
Evolutionary Phases 


For our single-star grids, we perform an additional 
post-processing step, which resamples the history out- 
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Figure 17. The state of compact objects is determined during the core-collapse step: if the mass of the pre-SN C/O core is 
less than 1.37 Mo, then we form a WD. Otherwise a SN will occur, and the state of the CO is determined based on its mass. 
For non-degenerate stars, the state is a combination of a surface-composition state, and a nuclear burning state. The former is 
decided solely by the presence or absence of hydrogen at the surface, whereas the burning state depends on whether a species 
(H, He or C) has been depleted from the core, and the burning of which species still contributes to the nuclear luminosity (Lnuc). 


put of the MESA code in a way that facilitates the inter- 
polation of entire evolutionary tracks. This is necessary 
for the computations described in Section 8.1. Using 
the method from Dotter (2016), we assign equivalent 
evolutionary phases (EEPs) throughout the evolution 
of a star. This method designates primary EEPs to ma- 
jor structural changes to a star (e.g., He ignition), and 
regularly spaced secondary EEPs in between. Primary 
EEPs are extracted directly from the computed stellar 
tracks. We then interpolate between the time-steps to 
identify every quantity of a star that we track at each 
secondary EEP. By applying this method to each of our 
single star tracks (both H-rich and He-rich) we can more 
easily interpolate within our single-star grids to find the 
quantities (e.g., radius, core mass) characterizing a star 
of any mass, and at any point throughout its evolution. 

The methodology described above is unfortunately not 
directly applicable to binary-star evolutionary tracks. 
For this interpolation method to work, the defined EEPs 
must be strictly ordered a priori. However, binary in- 
teractions, can happen at any point during the lifetime 
of a binary, often more than once, changing the order 
of EEPs in a non-tractable way. The interpolation of 


entire binary-star evolutionary tracks will be addressed 
in future releases of POSYDON. 


6.4. Downsampling of binary-star grids 


The evolutionary timesteps taken by MESA are typi- 
cally small, producing high-resolution binary histories 
and final profiles of the individual stars. To reduce the 
memory footprint of the data, and decrease computation 
times when modeling binary populations, we downsam- 
ple the binary tracks (i.e., keep a subset of the steps). 

For each individual run in a binary-star grid, we ob- 
tain from the MESA simulation the evolution of the bi- 
nary and individual stars’ parameters, as well as the 
post-processed quantities described in Section 6.2. For 
a total number of parameters m, the state of the binary 
is encoded in the m-dimensional vector h. The evolu- 
tion of the binary is a multivariate time-series given by 
h; = h(t;), with i = 1,--- , N, where N is the num- 
ber of steps, each corresponding to an age t;. Before the 
downsampling, the independent variable (age), and non- 
physical parameters (e.g., model number in MESA) are 
excluded from h, while all other parameters are rescaled 
linearly from 0 to 1. 
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Figure 18. The evolution of a 37.6 Mo star with a 22.6 Mo 
companion with an initial Po, of 6.82days. We compare 
the complete track provided by MESA (orange) comprised of 
3412 steps to our downsampled track (black dots) containing 
on 122 steps for the binary’s orbital period (top panel) and 
each stars’ radius (bottom two panels). In this particular 
case, the compression ratio is ~155, but the ratio varies from 
star to star and depends on which parameters are accounted 
for by the algorithm. The downsampling algorithm captures 
even the rapid variations seen in the two stars’ radii between 
4.8 Myr and 5.3 Myr (shown in the insets). 


The downsampling algorithm selects a subset of the 
original steps, so that if interpolated at the original 
timesteps t; the interpolation absolute scaled error is 
below a chosen threshold e: 


ei = ||h; — hy|| < €, (14) 


0.04 
All columns 
Selected columns 
0.03 
E 
e 0.02 
[au 


Compression ratio 


00.001 0.003 0.010 0.030 0.100 
Downsampling threshold 


Figure 19. The accuracy (Poo(r); top panel) and compres- 
sion ratio (bottom panel) of our downsampling algorithm as 
applied to the HMS-HMS binary grid, as a function of the 
downsampling threshold. We compare the performance of 
the algorithm when it is applied to all output columns in the 
data (circle markers) or only a selected list of columns (x; see 
excluded parameters annotated with an asterisk in Table 5 
and Table 6). As the downsampling threshold e increases 
from 107? to 107}, the compression ratio dramatically im- 
proves, but at the cost of the accuracy (Pos; represents the 
average of the 90-th percentile of the interpolation relative 
errors across all runs and either all or selected parameters). 
For our grids we use only selected columns with e = 107', 
which gives us a compression ratio of 726, while limiting 
any errors that could be accrued from this process to within 
a few per cent. 


where h; is the interpolated point using the age as the 
independent variable. We use linear interpolation, 
; ti — t 


hj = h; t : i (hj h;) , (15) 
titi — tj 


where j and j+1 are adjacent steps of the downsampled 
time-series so that tj < t; < tj44. 

The search for the steps in the downsampled data is 
performed as follows. Initially, we include only the first 
(hi) and last (hy) points. Then we search for the in- 
termediate point (h; where 1 < j < N) with maximum 
interpolation error. If this error is below the threshold 
(ej « €), then the algorithm has finished; otherwise it 
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includes this point and continues the search in the two 
parts of the time-series before (hi to h;) and after (h; 
to hy) the intermediate point. The process continues 
until all original points are well approximated by the in- 
terpolation of the selected subset of steps, i.e., e; « e, Vi. 

Additionally, we apply this downsampling method to 
the stellar profile data, following the exact procedure 
outlined above with the mass coordinate as the inde- 
pendent variable. A shell is kept not only when the in- 
terpolation error exceeds the predefined threshold, but 
also when the adjacent shells that are kept have differ- 
ence in mass is larger then 0.596 of the total stellar mass. 

'To demonstrate the validity of our method, Figure 18 
shows an example of the downsampling of a track using 
the same interpolation error threshold e — 0.1 we use 
for our grids. However, here we apply it to only two pa- 
rameters (orbital period and radius of secondary star) 
for visualization purposes; when the algorithm operates 
in a higher-dimensional space, it retains a large frac- 
tion of the initial points to capture the overall shape, 
making it hard to inspect its performance through two- 
dimensional plots. The downsampled version of the 
track is able to follow even the rapid oscillations oc- 
curring during the late-stage evolution of this particular 
binary. 

The choice of € is a balance between data compression 
and interpolation accuracy, a trade-off we demonstrate 
explicitly for our model grid composed of two H-rich 
stars in Figure 19. For our grids in POSYDON, we set the 
downsampling threshold to 0.1 and enforce it only for a 
list of 22 columns from the simulation output (see ex- 
cluded parameters annotated with an asterisk in Table 5, 
Table 6, and Table 7). This results to a compression fac- 
tor of ~26 with respect to original simulation data, but 
still sufficiently high accuracy with respect to the orig- 
inal grid. The final size of the three binary-star grids, 
after downsampling, is ~ 9.3 GB. 


7. OUR CLASSIFICATION AND INTERPOLATION 
APPROACH 


Even after our various stages of post-processing, we 
cannot use the grids of binary-star simulations within 
POSYDON as is for modeling populations (except if we fol- 
low a nearest-neighbor matching approach). While our 
binary-star simulations have only been run for a select, 
finite combination of initial masses and orbital periods, 
BPS requires us to have the capability to evolve a bi- 
nary anywhere within the domain of interest. We solve 
this problem in two steps. First, we apply a classifi- 
cation method to each of our grids to identify regions 
that undergo qualitatively different classes of evolution. 
Then we separately apply an interpolation method to 


each class to calculate stellar and binary properties. We 
describe the details of those methods in Section 7.2 and 
Section 7.3, respectively. 


7.1. Transformation and Rescaling of Grid Data 


For classification and interpolation purposes, we can 
interpret each of our binary grids as a data set that com- 
prises N input binaries, (x, Y ,, along with its cor- 
responding scalar {yn}; and class {z,,}_, targets. 
'These sets can be columnwise stacked into the matrices 
X € R?*N, Y e RM*N (where M is the number of out- 
put quantities), and Z of dimension 4 x N (because we 
classify each run into one of four broadly defined cate- 
gories; see Section 7.2). More specifically, each x, € R? 
contains the initial masses and orbital period of the n-th 
binary in the grid, whereas y,, and z,, denote the collec- 
tion of final binary and single-star quantities, and their 
associated classification, respectively. The N runs are 
distributed in a uniform three-dimensional mesh, either 
on a linear or logarithmic scale following the description 
provided in Section 5. This uniform grid constitutes an 
initial and naive way of thoroughly covering the parame- 
ter space, which is feasible due to the low dimensionality 
of the data. 

A convenient preprocessing of the data is crucial for 
both interpolation and classification. We apply a se- 
ries of non-linear and linear transformations to numeric 
data. Choosing the optimal transformation depends on 
the task (interpolation or classification), the method 
used for each task, and whether we are dealing with 
an input or an output quantity. 

First, we consider a non-linear transformation of the 
data using the logarithm: inputs can be transformed 
as logx; and targets as logy; or log(—y;) if y; < 0. 
Classification accuracy will improve when our algorithm 
uses the logarithm of the inputs for data sampled evenly 
in log-space. The effect on interpolation is different: 
e.g., a linear interpolation in the log-space results in 
a non-linear interpolation on the untransformed space. 
'This is similar to an approach where a non-linear space 
is transformed through a kernel to a space in which a 
linear model allows for modeling behavior appropriately. 

We automatically choose whether to apply a logarith- 
mic scaling using a cross-validation scheme. The optimal 
scaling for both inputs and outputs is chosen using the 
lowest relative error out of all the feasible scalings that 
could be applied to the given variable. 

As a second step, we apply a min-max scaling to the 
inputs so that the transformed features x are confined 
to the range [—1, 1], and we standardize the outputs such 
that they have zero-mean and unit variance, 

t yi — Yi 


(Yi Yi 1 
yi "m" (16) 
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Figure 20. The optimal k in our kNN classification scheme 
for each of our three grids. The bACC is calculated using 10- 
fold MC cross-validation, and averages the statistical recall 
of each of our four classes described in Section 7.2. The 
highest bACC for each grid is provided at the top of each 
panel and is indicated by the vertical lines; we use these ks 
when classifying our grids for running populations. Although 
we use the optimal k for each grid, our results are relatively 
insensitive to that choice. 


The choice of scaling for the inputs derives from the uni- 
form nature of the input grid data. Although it is pos- 
sible that we sampled our data in a non-optimal way, 
in practice we find the best results occur when our data 
scaling follows our grid sampling. In the case of the in- 
terpolated quantities, standarization produces improved 
metrics, particularly because it is less sensitive to out- 
liers. 


7.2. Classification of our Grids 


Accurate classification is a critical aspect of the 
POSYDON approach to evolving binary systems. There- 
fore, we separate our binaries into four categories based 
on their mass-transfer histories. The categories are: sta- 


ble mass transfer, unstable mass transfer, binaries that 
never interact, and those in RLO at ZAMS (Section 6.2). 
In addition to using their mass-transfer history, we could 
further segregate binaries into more refined classes; how- 
ever, we find this to be currently unnecessary, and we 
can accurately interpolate our binaries given these four 
broad classes. 

For each of our three binary-star grids, we generate a 
classification object that determines which of the four 
previously defined outcomes will be the result of a bi- 
nary with any particular combination of two masses and 
orbital period. In this first version of POSYDON we use 
a k-nearest neighbors (KNN) classifier, a simple and ro- 
bust classifier that achieves high precision in this task. 
We use the Euclidean distance as distance metric for the 
transformed input grid and weight each neighbor in the 
neighborhood proportionally to their inverse distance. 

We optimize the number of nearest neighbors we 
use by applying a Monte Carlo (MC) cross-validation 
scheme and selecting the k that produces a higher bal- 
anced accuracy (bACC). The bACC metric averages the 
statistical recall for each class (recall is the number of 
true positives divided by the combined number of true 
positives and false negatives) to produce a metric that 
accounts for any imbalances between classes. 

Figure 20 shows the average cross-validation perfor- 
mance for our three grids in terms of bACC as a function 
of the number of neighbors in the kNN classifier starting 
from k = 1 and highlights the location of the optimum. 
We train our final classifier on our regularly spaced grids 
using the optimal k, listed at the top of each panel in 
Figure 20 and indicated by a vertical black line. 

Figure 21 shows our classifier applied to one slice in 
each of our three grids, with the colors indicating dif- 
ferent regions. Overlaid grey contours, pronounced near 
class boundaries, indicate classification uncertainty. We 
ignore the no mass-transfer class for the grid of H-rich 
stars with CO companions, as this grid only applies to 
interacting binaries. 

To evaluate the accuracy of our classifiers, we use the 
validation data set associated with each of our three bi- 
nary grids. For each grid this validation set is comprised 
of binaries randomly sampled with the same range and 
scale, linear or logarithmic, as its training counterpart. 
Each of the three validation sets contain 3000 samples 
which roughly represent ~5-12% of the number of bi- 
naries in the regular grid. By applying our classifiers to 
the same initial values as those of our validation bina- 
ries, we can evaluate the accuracy of our classifiers. In 
Figure 21 our validation data is indicated by points (cor- 
rectly classified) and crosses (incorrectly classified). It 
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Figure 21. Decision boundaries of the kNN classifiers for a single slice in each of the three grids as a function of the primary 
star's mass on the horizontal axis and the orbital period on the vertical axis (the choice of q or Mco for each slice is indicated 
in each panel’s title). Shaded gray regions overlaid onto class regions represent the confidence of the classifier in that region. 
Points on top of the decision boundaries represent the validation data, where the edge color of each point shows the ground truth 
of the given point, and the fill color shows the classifier's prediction. Only in rare circumstances and only near classification 
boundaries does our classifier make incorrect predictions for our validation set. 
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Figure 22. Confusion matrices for each of our three binary grids. Each value at grid cell ci; represents the fraction of binaries 
that belong to class j (vertical axis) and were classified as class i (horizontal axis). Each row is normalized so that the sum of 
each row is 1, and the color of each cell indicates the magnitude of the value in the cell. Accuracies are all above 90%, with the 


exception of unstable mass transfer for the He-rich star with a CO companion. 


is evident that incorrectly classified validation binaries 
are very rare. 

To evaluate the quantitative accuracy of our classifier, 
we provide a confusion matrix for each of three grids 
in Figure 22. Diagonal squares indicate the fraction of 
systems that were correctly classified, while off-diagonal 
squares indicate the fraction of incorrectly classified sys- 


tems. The matrices are calculated such that each row 
sums to unity. All classes in all grids have an accu- 
racy in excess of 90%, often much more so, except for 
unstable mass transfer for binaries with a He-rich star 
and a CO companion. Examination of two slices of this 
grid in Figure 15 shows that the unstable mass-transfer 
class comprises a relatively small portion of the overall 
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grid, existing at small orbital periods, small CO masses 
and large companion masses. Reliable classification of 
small classes is difficult, but improving our classification 
accuracy will be a focus of future efforts (Section 11). 


7.3. Interpolation of our Grids 


Once classified based on their mass-transfer character- 
istics we separately interpolate binaries falling into each 
class for each of our three binary-star simulation grids. 
We only interpolate quantities for three of our binary 
classes, since those binaries overfilling their Roche lobe 
at ZA(He)MS are dismissed. 

We use an N-dimensional (where N is the number of 
binaries) linear interpolation: the data is divided into a 
set of N-simplices, tetrahedra in our three-dimensional 
data, by means of a Delaunay triangulation (which is 
not unique given the regular structure of our grids). The 
interpolated value for a given point corresponds to the 
value at the hyperplane that passes through the vertices 
of the simplex which contains the point. The choice 
of whether to apply a non-linear transformation on y;, 
logy; depends directly on that magnitude. For each 
output magnitude we select the optimal scaling via MC 
cross-validation with x iterations and p% of test data 
comparing the average relative error across iterations. 
'The final interpolator is trained using all binaries within 
a particular class for each grid. 

'The linear interpolation method is not capable of ex- 
trapolation: the value for any point which lies outside 
the convex hull defined by the constructed Delaunay 
triangulation will be undetermined. Although we are 
not, in general, interested in interpolating outside the 
training grid, there will be a small region between the 
convex hull of the linear interpolation and the decision 
boundary provided by the classifier where we still want 
to obtain system properties. For this small sliver of pa- 
rameter space, we adopt values of the nearest point in 
parameter space of the same class. This is a problematic 
region where the probability of belonging to the inter- 
polable class will be low, expressing the uncertainty we 
have about those binaries with the current resolution 
of the grids. We are currently exploring a method of 
tackling this problem by incorporating new simulations 
along the decision surface, identified using an active- 
learning scheme (Rocha et al. 2022). 


7.4. Ensuring Physical Congruity of Interpolated 
Values 


'The linear interpolation method described here treats 
each feature independently without preserving possible 
physical correlations. However, the interpolated results 
may produce incongruous quantities within a resultant 
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Figure 23. The convective envelope radius Reony.reg. and 
width Deonv.reg. need to simultaneously agree with three sep- 
arate inequalities defined in Table 4. The feasibility region 
(dark grey) represents the overlap of all three inequalities 
decomposed from the constraint. When our interpolation 
method proposed a point outside the feasibility region, we 
reassign it to a new point determined from the intersection 
of the border of the feasibility region and the line drawn 
from the centroid of the region to the constraint-violating, 
proposed point. 


star. For example, the Stefan—Boltzmann law connect- 
ing the luminosity, radius, and effective temperature of 
a star might not hold for an interpolated binary. As an- 
other example, the He-core mass of a star must always 
be less than the star’s total mass. We have carefully 
identified a number of physical constraints within the 
quantities that we are interpolating that must be satis- 
fied by any realistic star, each of which we list in Table 4. 

To address this issue, we process the interpolated 
quantities for a given binary so that they respect this 
list of constraints, which fall into one of three different 
Types depending on the basis of the corrective action re- 
quired. When quantities are connected via an equation 
(Type I), then the interpolated value of one parameter 
is ignored and inferred by solving the equation on the 
interpolated values of the remaining parameters. The 
Stefan—Boltzmann equation provides an example of a 
Type I constraint: we only interpolate each stars’ R and 
L, while Tg of the interpolated star is derived. In the 
case of inequalities between quantities (Type II), the 
quantity that must be less than another is limited by 
the value of the latter. Finally, there are cases where all 
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Constraint 


Relation 


Type I constraints: Equations 


Kepler’s Third Law 
Mass-transfer fraction? 
Stefan-Boltzmann Law 


Sum of nuclear luminosities 


a = [G (Mi + M3) P2,,/40?]*” 


r=1— Msys,2/ Mex 


Teg = (L/An B? ogg) ^ 


Imuc = Ly + Lye + Lz 


Type II constraints: Inequalities 


Mass-loss from the system from the vicinity of a star“? 
Core masses and radii 
Envelope masses and core radii? 


Mass, thickness and middle radius of the convective region 
(for tides) 


Remnant baryonic mass 


Moss < Mir and Msys,2 < Mir 
Mo/o-core < Mne-— core « M and Reo -core < Rue=core <R 
Mens « M and Reore<R 


Mconv.reg. « M and 
0 « Ficonv.xeg. — (1/2) Dcenv.reg. < Reonv.reg. + (1/2) Deonv.reg. < R 
Mrembar < M 


Type III constraints: constrained sum 


Central abundances 


Surface abundances 


Xeni 
Xs, H1 


XcN14 I X.c,016 ! Xc,other =1 
Xs, N14 + Xs,016 + Xs,other =1 


Xc,c12 
Xs,012 


Xc,He4 
Xs.He4 


'Table 4. Constraints that are ensured for our interpolated quantities. Type I relations are written such as the left-hand side 


indicates the quantity that is inferred from the rest. Notes: 


(2) This constraint is applied after constraint (b) since the mass-loss rate has to be less than the mass-transfer rate. Moreover, 


x is set to 1 if no mass is transferred. 


(9 There are four pairs of these quantities corresponding to the radii where the ! H fraction drops below 196, 10% and 30%, and 
finally where the ^He fraction drops below 1096 for pure He stars. 


quantities ought to add to a certain value (Type III). For 
instance, the fractional chemical abundances of a star's 
core must, by definition, sum to unity. We ensure these 
constraints are satisfied by normalizing our interpolated 
outputs. In one case described in the Table 4 footnotes, 
a parameter is subject to two separate constraints, in 
which case we are careful to apply them in the correct 
order. 

In the case of the constraint involving the interpo- 
lated quantities Reonv.reg. and Deony.reg., the middle 
point (= (Rs, conv.reg + Rb,conv.reg) /2) and the thick- 
ness (= Ri conv.reg — Rbp,conv.reg) Of the convective re- 
gion for the computation of tides, respectively, a spe- 
cial treatment is required. Both quantities must be 
positive and less than the star’s radius. However, con- 
straining them independently as in other Type II con- 
straints does not work as the inner and outer radius 
of the convective region must both be inside the star: 
0 « Ficonv.reg. m Deonv reg. /2 < Ficonv.reg. + Deonv.reg./2 < 
R. We decompose this relationship into three inequal- 
ities: Doonvreg:/2 > 0, Ficonv.reg. = D Shier 2; and 
Reonv.reg. + Deonv reg. /2 < R. In the Deonv.reg.~Reonv.reg. 
plane, the constraints form a feasibility region in the 
shape of a triangle with the vertices (0,0), (0, R), and 
(R, R/2), where R is a fixed values. If the constraints 
are violated, then the interpolated values lie outside of 
the triangle. The triangle’s centroid (R/3, R/2) and the 
point corresponding to the interpolated values define a 
line l. The intersection between / and the border of the 


triangle satisfies the constraint inequality, and is used to 
assign new values to the parameters. Figure 23 provides 
a pictorial representation of the algorithm. 

To assess how often our constraints defined in Table 4 
are violated in practice without imposing constraints, 
we interpolated 3000 binaries for each one of the three 
grids using the random initial conditions of their valida- 
tion sets (Section 7.5). For each binary, we checked all 
the constraints (two checks per binary system, and 23 
checks per non-degenerate companion star) and counted 
the violations. In the case of Type I constraints (equa- 
tions), we consider violation a relative error of more than 
0.001 in the inferred quantity. Overall, we found 57,035 
violations in the 274,560 checks (720.896) we performed. 
After applying the algorithm defined here, all violations 
were corrected. 


7.5. How Accurate are our Interpolation Methods? 


To assess the performance of the interpolation scheme, 
we use the same validation data sets that we used to 
evaluate our classification accuracy, described in Sec- 
tion 7.2. Our trained interpolators are applied to the 
same initial binary parameters as those of the three sets 
of 3000 binaries comprising our validation sets, one for 
each of our three binary grids. Since these binaries are 
not used in the training phase, the difference between 
this set and our interpolated predictions for them pro- 
vide an ideal comparison from which we can determine 
the accuracy of our methods. 
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Figure 24. Interpolation scheme accuracy for ten selected parameters when applied to our grid of two H-rich stars, as calculated 
using our set of validation binaries. We separate our sample by their different mass-transfer histories to independently evaluate 
their individual accuracy. Median relative errors (e) indicated by the horizontal lines in each distribution are typically 1% or 
lower for the stable mass-transfer and unstable mass-transfer cases (top panel) and the no mass-transfer case (green; bottom 
panel). Improving this accuracy will be a focus of future work. 
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Figure 25. The interpolation accuracy for the same 10 parameters as in Figure 24 for our grid of H-rich stars with a CO 
companion. Since we never use the models from this grid that do not undergo RLO, we do not evaluate the no mass-transfer 
binaries. In most cases, typical errors are 1% or better, but several of the distributions have tails extending towards larger e,. 


Figure 24 provides the accuracies for eleven selected 
binary and stellar parameters for our grid of two H- 
rich stars evolved from ZAMS. We have split these sam- 
ples by their mass-transfer histories so we can separately 
identify our algorithm’s accuracy for the stable mass 
transfer (red) and unstable mass transfer (blue) cases in 
the top panel and no mass transfer case (green) in the 
bottom panel. For nearly all parameters and all classes, 
our median errors are below 1%. Some parameters such 
as age and Jı are significantly more accurately interpo- 
lated, while others such as Mg/o-—core,1 may be some- 


what less accurate. The distributions are quite broad, 
suggesting that inaccuracies may exist when parame- 
ters show sharp variations as a function of input binary 
parameters. This is particularly apparent for our un- 
stable mass transfer channel, likely a result of the rel- 
atively smaller number of simulations that enter unsta- 
ble mass transfer and the varying evolutionary stages 
of the donor stars at the onset of the dynamical insta- 
bility. For instance, the relatively large error distribu- 
tion for Moo. core,1 for our unstable mass-transfer class 
is likely due to the rapid core growth during the giant 
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phase when donor stars typically enter dynamical insta- 
bility. Furthermore, despite their large error distribu- 
tions, some parameters, such as F4, have little impact 
on the evolution of a binary. Whether the primary star’s 
next evolutionary phase is a core collapse (in the case of 
the stable mass-transfer scenario) or a CE (in the case of 
unstable mass transfer) the mass at the outermost part 
of the star has little impact on the binary's outcome. 
Nevertheless, we plan to improve these accuracies with 
future enhancements to our interpolation schemes. 

In Figure 25 we provide analogous results for our sta- 
ble mass-transfer and unstable mass-transfer binaries for 
our grid of H-rich stars with a CO companion. Our mod- 
els tend to show larger variations in accuracy compared 
with our grid of two H-rich stars. Median errors tend 
to range from 196 to 1096 with certain parameters such 
as age and Mojo. core, performing noticeably worse. At 
the same time certain parameters like M» are very accu- 
rately determined, as these parameters vary during the 
evolution of the binaries in this grid (CO companions 
in this grid typically accrete little mass). One ought to 
consider the importance of each parameter when evalu- 
ating the accuracy of our models. For instance for stars 
going through unstable MT Mye is à much more im- 
portant parameter than Mc,o,; and Yurt has no impact 
on a binary's future evolution. 

Finally, Figure 26 shows the accuracies for our grid of 
He-rich stars with CO companions. Any edges in the 
distributions are genuine representations of the under- 
lying data. Our trained interpolators provide the most 
accurate predictions of our three grids; median accura- 
cies are typically between 0.196 and 196 for the stable 
mass transfer and unstable mass transfer classes, and 
somewhat better for the no mass transfer class. 

'The accuracies provided in Figures 24, 25, and 26 all 
refer to the data sets and associated interpolation ob- 
jects provided in v1.0 of POSYDON. One could use the 
POSYDON infrastructure to evolve larger numbers of bi- 
naries than we have provided along with v1.0, which 
would improve our interpolation accuracy. A focus on 
regions where our interpolation methods are least accu- 
rate would provide the largest benefit. Using a com- 
bination of active-learning techniques, more complex 
machine-learning algorithms, and much more computa- 
tion time, we expect that future versions of POSYDON will 
only exhibit substantially improved classification and in- 
terpolation accuracies (Section 11). 


7.6. Limitations of our Approach 


There are a few limitations of our approach. First, 
our approach first classifies the binary’s type and sub- 
sequently performs interpolation. The effect of such a 


technique is that by performing two optimization prob- 
lems, the second of which relies on the first, it is possible 
to propagate error throughout the pipeline. Treating the 
entire problem as one optimization problem has the po- 
tential to reduce error. 

Additionally, we transform the grid space by logarith- 
mic transformations before performing linear interpola- 
tion, which results in a non-linear interpolation model. 
Such an approach is similar to using a kernel, where a 
space is transformed through a kernel function to a space 
in which a linear model allows for accurate modeling of 
the behavior of the space (Theodoridis & Koutroumbas 
2009). A more systematic approach is to consider a ker- 
nelized interpolation approach, by applying kernel selec- 
tion techniques. In the case of a Gaussian process, for 
example, we may consider a whole family of functions 
which are specified by a kernel function, to allow for 
more flexibility on the prior belief of the space (MacKay 
2003). 

Finally, upon finding an interpolated value, we physi- 
cally enforce the constraints, as detailed in Section 7.4. 
However, in principle such a technique, which does not 
consider the constraints in the optimization objective it- 
self does not guarantee an optimal solution subject to 
the constraints. One way to incorporate the constraints 
in our model is to add a regularizer term in our loss func- 
tion to enforce the constraints (Ivezić et al. 2014), i.e., 
the loss function balloons when constraints are violated. 


8. EVOLUTIONARY PROCESSES SEPARATE 
FROM SINGLE- AND BINARY-STAR MODEL 
GRIDS 


Besides computing, processing, classifying, and inter- 
polating the five separate grids of single- and binary-star 
models, additional steps are required to follow the com- 
plete evolution of a stellar binary from ZAMS to dou- 
ble CO formation (and potentially its merger). These 
are defined by three separate processes: orbital evolu- 
tion in eccentric, detached binaries, CE evolution, and 
stellar core-collapse. While the latter two are standard 
elements of BPS codes, the need for the former requires 
some explanation. Binaries are intrinsically eccentric af- 
ter a SN occurs, yet our pre-calculated grids of binary- 
star models are initiated with circular orbits. Including 
eccentricity as an input to our MESA models would add an 
additional dimension to our simulation grids, challeng- 
ing our computational capabilities. Furthermore, self- 
consistently modeling binary mass transfer along with 
stellar evolution and tides in eccentric orbits is an ac- 
tive area of study (Sepinsky et al. 2007, 2009, 2010; 
Dosopoulou & Kalogera 2016a,b; Hamers & Dosopoulou 
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Figure 26. The interpolation accuracy for the same 10 parameters as in Figure 24 and Figure 25, but for our grid of He- 
stars with CO companions. As with our other two grids, this grid has typical median errors below 1%; however, tails of the 
distribution extend towards larger e;, especially for the stable mass-transfer and unstable mass-transfer cases (top panel). Our 
no mass-transfer case (bottom panel) is much more accurate than the corresponding binaries in our grid of two H-rich stars in 
Figure 24. The apparent truncations in the distributions (e.g., Mc/o-core,1 and log,9(J1)) are genuine representations of the 


data. 


2019), and to date no detailed binary evolution grids 
have included initially non-circular binaries. 

Nevertheless, in a detached binary, tidal forces cause 
an eccentric binary to both circularize and synchronize, 
an effect that must be taken into account, along with 
other orbital angular-momentum loss processes (e.g., 
wind mass loss, gravitational radiation and magnetic 
breaking). To specifically address this, we evolve bi- 
naries after the a SN event using a separate process de- 
scribed in Section 8.1. We only switch back to using 
the pre-calculated grid of binary-star models once RLO 
occurs. 

In the current version of POSYDON, binaries that suc- 
cessfully exit from a CE phase initiated by two non- 
degenerate stars, are also modeled following the process 
described in Section 8.1 (see also Figure 1). These are bi- 
naries consisted of a H-rich and a He-rich, or two He-rich 
stars in a close circular orbit. We follow their evolution 
as detached binaries until one of the two stars reaches 
core-collapse. In a small fraction of post-CE binaries 
(~ 0.3% of the total population) typically consisted of 
a low-mass (S 4Mo) He-rich and a H-rich MS star, 
the He-rich star overflows its Roche-lobe as it expands 
to become a giant He star, and initiates mass transfer 
onto the H-rich MS star. Since we do not have a grid 


of detailed binary-star models covering this part of the 
parameter space, we cannot follow further the evolution 
of these binaries. We plan to address this in future ver- 
sions of POSYDON. Finally, the evolution of all binaries 
that are in orbits wider than what is covered by our 
grids of detailed binary-star models, and thus will never 
initiate RLO, is also followed as described in Section 8.1. 

In the subsequent sections, we provide details about 
how we evolve binaries through an eccentric, detached 
phase, a CE evolution phase, and core collapse. 


8.1. Evolution of eccentric, detached binaries 


8.1.1. Matching with a single-star track 


Even though a non-degenerate star in a detached bi- 
nary is influenced by its companion (for instance, due 
to tides), we are making the assumption that so long 
as RLO is not occurring, our single-star, non-rotating 
models provide reasonable approximations for the evo- 
lution of non-degenerate stars in these detached binaries. 
We first match the non-degenerate star in the detached 
binary system with the closest model (searching across 
different masses at all ages) from our single-star (both 
H-rich and He-rich) evolutionary tracks. The matching 
is achieved by minimizing the sum of the squares of key 
parameters describing the structure of a star. These pa- 
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Figure 27. Distributions of the difference between the matching point of a single star model in the beginning of the detached 
step and the values from the previous step, for various physical quantities of a non-degenerate star. We show the relative 
difference in mass AM, He-core mass AMye~—core and logarithm of the radius A log,o(R/ Ro), as well as the difference in the 
He central and surface abundances Yeenter and Ysurf, respectively. Vertical, dashed lines from left to right delineate the 5-th and 
95-th percentiles of the distributions. The distributions all suggest that the non-degenerate stars in detached binaries can be 


accurately matched to a single star model. 


rameters differ depending on the evolutionary phase of 
the star. 

For the matching process, we distinguish among: (i) 
MS stars that still have H in their core (with central H 
mass abundance Xcenter > 0.01), (ii) stars that evolved 
off the MS and retain even a thin H-rich envelope (post- 
MS, with Xcenter < 0.01 and surface H mass abundance 
Xsurt > 0.01), and (iii) evolved stripped stars that are 
effectively a H-deficient core (Xsurf < 0.01). In case (i) 
the parameters whose differences are minimized are the 
total mass of the star, its central H abundance X center 
and the radius of the MS star. After the first core col- 
lapse of the system occurs, most binary companions are 
in this state, being the initially less massive secondaries 
that evolved slower than the primary that has formed a 
CO. For case (ii), we replace the central H abundance 
with the He one (Yeenter) and the radius with the mass 
of the fully developed He core. In case (iii), for stripped 
He-rich stars, we use as minimized parameters the He- 
core mass of the star (which is equal to its total mass), 
the radius, and its center He mass abundance. 

We normalize the chosen quantities, such that they 
have similar weighting in the minimization process. The 
normalization factors are chosen from typical ranges 
of each parameter for the stars that we focus on: 
20 Mofor the total mass of MS or post-MS stars; 10 Mo 
for He-core masses of stripped He-rich stars, 2.0 for 
log;59(R/ Ro); and 1.0 for chemical mass-fraction abun- 
dances. 

We quantify the quality of the matching by calculat- 
ing the difference of various quantities from the pre- 
vious step. In Figure 27, we see that the difference 
between the new interpolated total mass in the be- 


ginning of the detached step and its previous value, 
AM = Myatch — Mprev,step, is typically better than 
0.2 Mo. He-core masses are matched even more pre- 
cisely, to within 0.05 Mc. Other parameters (we show 
R, Yeenter; Ysurt in Figure 27) are also closely matched, 
justifying our assumption that the non-degenerate star 
in a detached binary can be accurately represented by a 
single star model, at least so long as it remains detached. 

A CO component in the binary system is treated as a 
point mass and does not need a matching process. We 
also keep the orbital separation and eccentricity con- 
stant in this transition, and thus any small difference 
in the matched star’s mass from the previous step re- 
sults in a small relative change in the orbital period. 
In addition, by conserving the spin angular momentum 
on non-degenerate stars during the matching step, we 
can determine their initial angular frequency 2 at the 
beginning of the detached step. 


8.1.2. Further evolution of an eccentric detached binary 
system 


Once matched with single star models, we evolve the 
stars in detached binaries as essentially single stars, ac- 
counting for their effects on the binary's orbit. For a 
non-degenerate star, its parameters (e.g., mass, radius, 
moment of inertia) are evolved according to its interpo- 
lated stellar track. At the same time, its spin 2 as well 
as the system's a and e are evolved solving a set of cou- 
pled ordinary differential equations that describe their 
rate of change due to wind mass loss, tides, magnetic 
braking, and gravitational radiation. That we assume 
the back reaction of each of these effects does not sig- 
nificantly impact the internal structure of either star so 
that our single-star models are sufficiently accurate. For 
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instance, although we follow each star’s spin using the 
moment of inertia of the single-star models, we cannot 
account during this phase for the star’s internal differ- 
ential rotation and effects such as rotational mixing. 

This approach can only handle scenarios where no 
RLO mass transfer takes place between the two stars; as 
soon as a H-rich star enters RLO, we stop the binary’s 
evolution and transition to our grid of MESA mass trans- 
fer simulations described in Section 5.6. Binaries are 
assumed to circularize instantaneously upon RLO with 
an orbital separation equal to the binary’s separation 
at periastron. Alternatively, we also allow for a user to 
choose to circularize the orbit assuming angular momen- 
tum is conserved.? Likewise, this step of evolution also 
ends if a non-degenerate star reaches the end of its life, 
in which case the binary is sent to a step that handles 
core collapse (Section 8.3). 

A third stopping condition exists for binaries consisted 
of two COs, which merge due to gravitational wave ra- 
diation. Note that when modeling two COs, only ef- 
fects due to gravitational wave radiation contribute to 
orbital evolution. We calculate the orbital decay until 
the merger or the maximum simulation time is reached. 

Orbital evolution during the detached step is due to a 
combination of the relevant pieces of physics, which we 
assume have additive effect: 


a= wind T tides,1 sr tides,2 =F agar; 

é- Étides,1 T Étides,2 s EGR; 

Qi — Ouwind,1 =F QT ae Ouaes,1 + Ops 
OQ = Owind,2 T Öintertia,2 HF Crides, 2 + Don 


The orbital separation, eccentricity and stellar spins are 
evolved using a set of self-consistent, coupled equations. 
We describe each of the terms in Eq. (17)-(20) below. 

Mass Loss: We ignore mass accretion onto a star 
(either non-degenerate or CO) from a companion star's 
wind. So non-degenerate stars will only lose mass due 
to their own stellar winds, with the mass lost carrying 
away the specific orbital angular momentum of the mass- 
loosing star (Jeans-mode mass loss; for a review, see 
Tauris & van den Heuvel 2006): 


w,1 T w,2 
L, 21 
Mı + M5 M 


awind = —@ 


For CO binary components, My = 0, while in general 
for spherical, isotropic fast winds, the orbit-averaged e 


6 This option results to circularized orbits where the star does not 
fill its Roche lobe anymore. In this case, we allow for the start 
to further evolve until it fills its Roche lobe again, but without 
changing the orbit. 


due to winds is zero. We discuss the effect of mass loss 
on stellar spin later in this section. 

Tides: Changes in the orbit's period and eccentricity 
due to tidal forces are described by a set of ordinary 
differential equations, according to Hut (1981). In order 
to be able to compute tidal spin-orbit coupling, we treat 
the donor star of mass M, radius R and of moment of 
inertia I, as a solid body rotating with angular velocity 
Q. The initial angular-momentum budget of the non- 
degenerate star is assumed to be the same as from the 
end of the previous step. 

'The change of the orbital separation due to tidal forces 
on the first star (subscript 1) is given by: 


kN. Mı(Mı + M35) / F3? a 
Atides,1 = 6 
1 (1 


T M? a — e2)5/2 


x |^ (2) - 1 - e? (e) 2 | (o (2) 
orb 

where Qorb = 2m-/P4, is the mean orbital angular 
velocity. "When both stars are non-degenerate, they 
each have their own contribution to the orbit's evolu- 
tion. Therefore an analogous equation exists providing 
Gtides,2, Where R; is replaced with Rə, the k/T term is 
calculated for the secondary star, and Mı and Ms» are 
switched. 

The k/T term in Eq. (22) depends on a star’s structure 
and the associated physical process of tidal dissipation. 
We calculated them separately for dynamical and equi- 
librium tides, in the same way as in our detailed, binary- 
star model grids (Section 4.1) described in Eq. (4) and 
Eq. (6), respectively. We apply the maximum of these 
two at each timestep. 

'The change of the orbital eccentricity and the stellar 
spin from tidal forces is also calculated as 


anc o car ( E) Mi +M) (R j e 
tides,1 = T i M2 a (1 e2)18/2 
11 3/2 Q 
x | (e?) - TA (1 — e) / fa (e) | "c 
Sue (2), G2) Co) aes 
tides,1 — T A Mı qi a (1— e?)° 


x Ẹ (8) - (1— &*? & (e) a . (uU) 


As in the Eq. (22) for àjaes, when the companion star is 
non-degenerate, €tides,2 and Okides,2 terms exist, which 
can be calculated by switching subscripts 1 and 2 in 
Eq. (23) and (24). The f;(e?), i = 1—5 terms in Eq. (22), 
(23), and (24) can all be found in Hut (1981). 
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Stellar Evolution: During the detached orbital evo- 
lution, we also take into account the change of stellar 
spin due to the evolution of the stars themselves. This 
includes spin down because of wind mass loss that car- 
ries away the specific angular momentum of the star’s 
surface, as well as changes in its spin due to the evolution 
of its moment of inertia due to changes in its internal 
structure, 

: 2 R2 0, ,. 0; 


Quin inertia,] — M I, 25 
d+ tia,1 3 L 1 Ñ 1 ( ) 


where I the rate of change of its moment of inertia. 
For binaries in which both stars are non-degenerate, an 
equation equivalent to Eq. (25) exists for Oandiinertia 2: 

Non-degenerate stars tend to spin down, due to their 
expansion and their wind mass loss. However, they may 
also be spun up in phases where they contract. There- 
fore, for numerical-stability reasons we artificially limit 
the second term in Eq. (25) to +100radyr~? (usually 
reached during a sudden contraction to form a WD). 
Although we take into account the effect of stellar spin 
on the orbit via tidal spin-orbit coupling, we do not 
include effects of spin on the stellar structure, such as 
stellar deformation, rotational mixing or rotationally- 
enhanced winds. 

Magnetic Braking: In case the binary contains low- 
mass non-degenerate stars, spin-down due to magnetic 
braking can become important i.e., the loss of spin an- 
gular momentum due to ionized material ejected from 
the star that is trapped in its own radial magnetic field. 
The spin-down rate is given by 

hoc 
where Tmb,1 is the torque calculated as in Eq. (36) of 
Rappaport et al. (1983), 


up, = (26) 


M R Ymb Q, /2 3 
Tmb,1 = —6.82x10°4dyncm ( -) ( -) ( 1/ z) , 


Mo Ro 1/day 

(27) 
with Ymb = 4 (Verbunt & Zwaan 1981). We apply the 
full torque to all non-degenerate stars below 1.3 Mc and 
assume no magnetic braking for stars above 1.5 Mo, 
with linear interpolation in-between. Again, for bina- 
ries with two low-mass, non-degenerate stars, an equa- 
tion equivalent to Eq. (26) exists for se 

Gravitational Wave Radiation: Finally, we take 
into account orbital changes due to gravitational radia- 
tion (Peters 1964; Junker & Schaefer 1992), 


2 vc on Lx 
OGR = 


15 (1 " e2)9/2 
1 G(M; + M3) 
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SR 15 G(M, + M3) ac? (1 — e2)? 
1 G(M; + M3) 
al- Oa ()], — qu 


where v = Mı M3/(Mij 4- M3)?, and the functions g;(€?), 
i — 1 — 4, are defined as: 


gi (e?) = (96 + 292e? + 37e*) (1 — e°), (30) 
g2 (e?) = (14008 + 4704v) + (80124 + 21560v) e? + 
(17325 + 10458) e* — 


1 
5 (5501 — 1036») ef, (31) 

gs (e?) = (304 + 121?) (1 — e°), (32) 

ga (e°) =8 (16705 + 4676v) + 12 (9082 + 2807v) e? — 
(25211 + 3388v) e^ . (33) 


These general-relativistic terms of orbital evolution 
are usually negligible apart from cases of close binaries. 
For binaries consisted of two COs, only Eq. (28)-(33) 
govern the evolution of the binary's orbit. 


8.2. Common-envelope evolution 


Binary interactions can lead to a dynamically unsta- 
ble mass-transfer phase (Ivanova et al. 2013, 2020). We 
have described in Section 4.2.4 all the conditions that are 
assumed to trigger an unstable mass-transfer episode: a 
maximum mass-transfer rate of 0.1 Mọ yr-!, Lz over- 
flow, a contact phase with a post-MS star, or exceeding 
the threshold of the trapping radius during accretion 
onto a CO. 

If the donor star that triggered the unstable phase 
is in its MS phase or is a stripped He star during its 
He core-burning phase, we assume that the two stars 
promptly merge, as no distinct core has formed yet in 
its interior. In v1.0 of POSYDON we do not follow the 
further evolution of stellar merger products. 

For all the other donor stellar states, a trigger of un- 
stable mass transfer is assumed to lead to a CE phase. If 
the donor has a H-rich envelope at the beginning of the 
phase, this envelope is considered to form the CE, inside 
of which the donor’s He-rich core and its binary com- 
panion will spiral-in. For stripped donors, the He-rich 
envelope engulfs the companion which spirals in around 
the donor's C/O core. In case the companion star also 
has a giant-like structure with a distinct core-envelope 
separation, (i.e., anything but a MS star, a He star in its 
He-MS, or a CO), then its envelope also may contribute 
in a (double) CE. 
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Figure 28. Evolution of Ace parameter of POSYDON single- 
star models of different initial masses. For these calculations, 
the assumed core-envelope boundary is located at the point 
where the H mass fraction drops below 1096. The triangle, 
diamond, and star markers represent the start of shell H 
burning, the start of core He burning, and the end of core 
He burning, respectively. 


The outcome of the CE phase is calculated using the 
QcE-AcE prescription (Webbink 1984; Livio & Soker 
1988), which equates a fraction agg of the orbital energy 
released during the spiral-in with the binding energy of 
the CE. The acg parameter is set equal to 1 in the ex- 
ample population runs shown in Section 10, following 
previous population synthesis works (e.g., Hurley et al. 
2002), but is, in general, a free parameter in POSYDON. 

The parameter Ace has been introduced to 
parametrize the binding energy of the envelope us- 
ing the total stellar radius and mass (de Kool 1990). 
In POSYDON, Aog values are calculated from the de- 
tailed stellar profile of the donor star at the beginning 
of CE (or of both stars, in case of a double CE). This 
is an important quantitative improvement of POSYDON, 
compared to pBPS codes. The latter need to adopt 
Acrg-value fits from the literature, based on single-star 
models with often inconsistent stellar-physics, and ap- 
ply them to post-interacting stars. In our common CE 
energy calculation, we integrate both the gravitational 
and the internal energy of the envelope from the de- 
tailed stellar structure profiles of our binary models, 
subtracting the recombination energy from the internal 
energy. In Figure 28 we show the Ace for a few example 
single-star POSYDON models. The parameter Acg tends 
to decrease as the star evolves and expands as a giant. 
However, for initially very massive (>, 20 Mo) stars that 
strip their H-rich layers due to their own wind mass loss, 
Ace increases again. We find comparable trends and 


values with other works that study the detailed stel- 
lar structures of giant stars (e.g., Kruckow et al. 2016; 
Klencki et al. 2021). 

The binding energy of the envelope (and thus the out- 
come of the CE phase) is sensitive to the exact assumed 
core-envelope boundary of the donor, as the deeper en- 
velope layers tend to be the most tightly bound (Dewi & 
Tauris 2000; Ivanova et al. 2013; Fragos et al. 2019). For 
this reason, we allow for different core-envelope bound- 
aries, defined for H-rich stars as the outermost layer 
where the H mass fraction drops below 0.3, 0.1 (default), 
and 0.01 and for stripped-stars when the sum of H and 
He drops below 0.1. 

Given the properties of the binary at the onset of the 
CE, the assumed acg value and the estimated Acg, one 
can calculate how much a binary's orbit shrinks in order 
for the released orbital energy to fully unbind the CE. 
The final post-CE orbital separation aposc,cg is given by 
solving (Webbink 1984): 


G Maon,core Macc 


Zápost,CE 


G Maon Macc un 


2apre,CE 


G Maon Maon,env 

QGCEACE aon 

(34) 
where Maon, Maon,core, and Maon,env are the total, core 
and envelope masses of the donor star, Rao, the donor 
star's radius, Macc the mass of the accreting star, and 
Gpre,CE the orbital separation at the onset of the CE. 
If the final estimated apost,ce is such that neither the 
accreting star nor the stripped core of the primary star 
are filling their respective Roche lobes, then the CE is 
consider successful and results in a detached, circular, 
tight binary. Alternatively, the binary is assumed to 
merge, and its evolution is not further followed in v1.0 
of POSYDON. 

One complication with the flexibility we offer regard- 
ing the core-envelope boundary definition is that the 
post-CE stripped donor star might still contain some H 
in its outer layers, while in the next evolutionary steps 
we assume that the H envelope is fully removed. Ex- 
actly how much H remains depends on a user's choice 
of 0.01, 0.1, or 0.3 for a fractional H abundance when 
defining the core-mass boundary. We account for this 
inconsistency by assuming that either the remaining H- 
rich layers are either removed by stellar winds or these 
layers re-expand after the CE and are removed via stable 
mass transfer (e.g., Fragos et al. 2019). Both assump- 
tions result in slight corrections to the post-CE donor 
masses and orbital separation. Although the former as- 
sumption is the default one in POSYDON, we find they 
both lead to correction at the level of only a few per- 
cent, and thus the choice between the two is in practice 
inconsequential. 
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8.3. Core-collapse and compact-object formation 


The end fate of stars primarily depends on their 
masses. The most massive stars undergo all nuclear 
burning phases (hydrogen, helium, carbon, neon, oxy- 
gen, silicon) up to the formation of an iron core . The 
iron core keeps growing by silicon shell burning to a mass 
of around the Chandrasekhar mass limit ~ 1.44 Mo 
when electron degeneracy pressure can no longer sta- 
bilize the core and it collapses. This runaway process 
can lead to the explosion of a star in a SN or to a direct 
collapse into a BH and it is known as core-collapse SN 
(CCSN) (see Janka et al. 2007, for a review). 

Lower mass stars do not complete all nuclear burning 
phases. For stars which do not ignite oxygen but for 
which their He-cores masses are 1.4 Mo S Mue-core S 
2.5 Mo (Podsiadlowski et al. 2004) we assume a star 
collapses into a NS in an electron-capture SN (ECSN). 
POSYDON alternatively includes the option to determine 
whether a star undergoes an ECSN based on its C/O 
core mass: 1.37 Mo S; Mojo-cee S; 1.43 Mo (Tauris 
et al. 2015). Stars with core masses below the lower 
limit for ECSN evolve into white dwarfs. 

In Figure 29, we show, for a slice at a fixed initial mass- 
ratio q — 0.7 of the binary-star model grid composed of 
two H-rich stars, the core-collapse type as a function of 
initial orbital period primary star mass. The transition 
region between the ECSN and CCSN, occurs at ZAMS 
masses of ~ 8 Mọ (consistent with previous studies, e.g., 
Nomoto 1984; Jones et al. 2014), but depends somewhat 
on the initial Pop. 


8.3.1. Pulsational pair-instability SN 


During the post-carbon burning phase of massive stars 
(not modelled here), photons produced in the core can 
be energetic enough to produce electron-positron pairs, 
softening the equation of state and diminishing the pres- 
sure support of the core (Woosley et al. 2007, and refer- 
ences therein). In such stars the core rapidly contracts 
and the temperature increases, leading to explosive oxy- 
gen burning (e.g., Woosley & Heger 2015) that creates 
a series of energetic pulses which eject material from 
the star surface. This phenomenon of material ejection 
due to pulses is known as pulsational pair instability SN 
(PPISN) and occurs for stars with He-core masses in 
the range ~ [32,64] Mo, (Yoshida et al. 2016; Woosley 
2017; Marchant et al. 2019; Renzo et al. 2020). For more 
massive stars with He-core mass in ~ [61,124] Mo, the 
first pulse is so energetic that can unbind and destroy 
the whole star in a so-called pair instability SN (PISN; 
Fowler & Hoyle 1964; Rakavy & Shaviv 1967; Barkat 
et al. 1967), leaving no remnant behind. 
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Figure 29. The core-collapse type for the q = 0.7 slice 
of our binary-star model grid composed of two H-rich stars. 
We distinguish between WD formation, electron caption SN 
(ECSN) following Podsiadlowski et al. (2005), core-collapse 
SN (CCSN) and pair-pulsational instability SN (PPISN) fol- 
lowing Marchant et al. (2019). Models that did not reach 
the end of stellar evolution are indicated in black. 


To identify systems that will undergo PPISN and 
PISN, we adopt a polynomial fit, as implemented in 
Breivik et al. (2020, see Eq. 4 therein), to MESA single- 
star simulations (at Z = 0.1Z5) by Marchant et al. 
(2019, see Table 1 therein). This fitting formula is used 
to map the He core mass at carbon depletion in the 
range 31.99 Mo < Muecore < 61.10 Mo to the stellar 
mass collapsing to form the CO. We use PPISN models 
computed at 1/10-th solar metallicity, as it was shown 
that such a limit is independent of metallicity (Farmer 
et al. 2019), while highly dependent on the uncertain 
12C(a, )190 reaction rates (Farmer et al. 2020). In our 
case, these reaction rates follow Cyburt et al. (2010), 
consistent with the rates used by Marchant et al. (2019). 

In Figure 29, we can identify two systems that enter 
the regime of PPISN as they possess a He-core with 
mass slightly larger than 31.99 Mo at carbon depletion. 
Other mass ratio slices show a few more similar systems 
but without becoming statistically relevant. We expect 
PPISNe to be more present at sub-solar metallicities, 
as stellar wind-mass loss at Zo prevent the stars from 


THE POSYDON BINARY POPULATION SYNTHESIS CODE 41 


reaching carbon depletion with a He-core mass in the 
relevant mass range for PPISNe. 


8.3.2. Remnant baryonic mass 


In this version of POSYDON, we calculate the mass left 
behind by the collapse using different models: (i) direct 
collapse where all the stellar mass is conserved; (ii) fits 
to the results of two-dimensional core-collapse models 
of Fryer et al. (2012); (iii) nearest neighbor interpola- 
tions of the results of the detailed one-dimensional core- 
collapse models of Sukhbold et al. (2016), or (iv) with 
the explodability criteria of Patton & Sukhbold (2020). 
The last is our default option. 

Fryer et al. (2012) presents two mechanisms which 
are known as rapid and delayed based on how quickly 
convective instabilities are expected to grow after core 
bounce. The rapid prescription produces a mass gap 
between BHs and NSs by assuming strong convection 
which allows instabilities to grow quickly after core 
bounce, producing a more energetic SN explosion. In 
contrast, the delayed mechanism produces a continu- 
ous spectrum of compact remnant masses. Both pre- 
scriptions determine the baryonic mass of the com- 
pact remnant Mrembar given the pre-SN C/O core 
mass, Mco-—core- More precisely, Mc/o-—core determines 
whether the star explodes into a SN, and what fraction, 
ft», of the ejected mass falls back onto the CO. In the 
case that the star directly collapses to form a BH, fp = 
1. For the rapid prescription, direct collapse occurs for 
Mc/o-core 2 7.6 Mo while for the delayed prescriptions 
direct collapse occurs for Mcjo—core 2 11 Mo. 

In Sukhbold et al. (2016), the outcome of the col- 
lapse of their pre-SN models has been calibrated against 
the well-studied SN 1987A progenitor. We have imple- 
mented several of their SN engine calibrations, namely 
N20, $19.8, W15, W20 and W18, although for the train- 
ing of the initial-final interpolation we only consider the 
default N20 option, which is the most optimistic op- 
tion for successful explosions. In contrast to Fryer et al. 
(2012) results, Sukhbold et al. (2016) finds sharply vary- 
ing behavior between the initial star mass and the fi- 
nal core properties, linked to convective carbon-burning 
episodes occurring in the later evolutionary phases. This 
results in a region of the parameter space where the 
outcome of the collapse, i.e., NS and BH formation, ap- 
pears stochastic in its nature. To mitigate interpolation 
errors, we determine the remnant baryonic mass of a col- 
lapsing star using a nearest neighbor technique on the 
He-core mass at carbon depletion to map our stars to 
the Sukhbold et al. (2016) simulation results. 

In the Patton & Sukhbold (2020) prescription, the 
C/O core mass and the average carbon abundance of the 


core at carbon ignition are used to determine the explod- 
ability of the pre-SN core. For every single- and binary- 
star model in our grids, we store these two values, and 
by applying a k-nearest neighbour interpolation, with 
k — 5, we map to the explodability parameters M4 and 
L4 from Ertl et al. (2016), as described in Patton & 
Sukhbold (2020). These two explodability parameters 
allow us to infer whether a SN is successful, and, if so, 
we estimate the resulting NS mass to be approximately 
equal to M4. We assume that BHs are produced only 
from failed explosions which result in a direct collapse. 
Finally, for the Patton & Sukhbold (2020) prescription, 
we have implemented the same SN engine options as for 
(Sukhbold et al. 2016) with N20 as our default option, 
using the updated calibration from Ertl et al. (2020). 

In Figure 30, we show a comparison between the final 
CO state for the same grid slice as Figure 29, as pre- 
dicted by Fryer et al. (2012) delayed prescription com- 
pared and the Patton & Sukhbold (2020) prescription, 
based on the N20 engine. The differences between our 
choice of SN prescription are slight, but noticable when 
focusing on the NS/BH boundary. The Fryer et al. 
(2012) delayed prescription produces BHs for somewhat 
less massive stars, while the Patton & Sukhbold (2020) 
prescription shows a more variable boundary between 
NSs and BHs. 

In both Sukhbold et al. (2016) and Patton & Sukhbold 
(2020) prescriptions we assume fallback fractions of 
fm = 1 for BHs and fm = 0 for NSs. For Fryer et al. 
(2012) prescriptions the fallback fractions are computed 
explicitly, with the exception of NS ECSN where we as- 
sume fm = 0. 


8.3.3. CO gravitational mass 


'To convert the remnant baryonic mass to gravitational 
mass, we use the prescription by Zevin et al. (2020), 
which is an updated version of the one by Lattimer 
& Yahil (1989) based on the neutrino observations of 
SN 19874. This new conversion caps the maximum neu- 
trino mass loss to 0.5 Mo (C. Fryer, private communi- 
cation) and removes any artificial discontinuity in the 
mass spectrum between NS and BH formation (in the 
case of direct collapse or the Fryer et al. (2012) delayed 
mechanism) as 


20 Myembar 
1+0.3 
3 ( Mo 


1) Mo, 


Miembar = 0.5Mo otherwise. 
(35) 
If Maa, < 2.5 Mo we classify the CO as a NS other- 


wise as a BH. There is a large uncertainty in the ex- 


Miembar = Meray « 0.5 Mo 
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Figure 30. The CO state for the q = 0.7 grid slice of HMS-HMS MESA simulations in the initial primary mass-orbital period 
plane. We distinguish between the COs: white dwarf (WD), neutron star (NS) and black hole (BH) according to the legend. 
Models that did not reach the end of the stellar evolution are indicated in black. The two panels compare the Fryer et al. (2012) 


delayed core-collapse mechanism (left) with the outcome of Patton & Sukhbold (2020) N20 core-collapse engine (right). 


star which forms a proto-BH of mass 2.5 Mo. The mass 
lost in neutrinos during the formation of this proto-BH, 


act maximum NS mass and this range spans 2.0-2.7 Mo 
(Lattimer & Prakash 2010; Margalit & Metzger 2017; 


also carry away specific angular momentum equal of that 


Rezzolla et al. 2018; Ai et al. 2020; Shao et al. 2020; 
Lim et al. 2021; Miller et al. 2021; Raaijmakers et al. 


2021). In POSYDON, the maximum neutron star mass is 


collapsing central part of the core that forms the proto- 


BH. If Mrembar < Mstar we assume the ejected mass 


takes away the outer layers of the star. 


set to Myg* = 2.5 Mo (see discussion in Abbott et al. 


2020b, and references therein). 


'The angular momentum content of the infalling mate- 
rial can in principle support the formation of an accre- 


In the left panel of Figure 31, we show the grav- 
itational mass of the CO as predicted by Patton & 


Sukhbold (2020) N20 prescription for the same grid slice 


as Figure 30. 


tion disk. We consider a collapsing star to be a collection 


of shells with radius r, mass Mshey and angular velocity 


shell that falls one by one onto the central BH. A shell 


of mass is accreted by the BH once it reaches the BH's 
event horizon. The specific angular momentum of the 


Q 


Finally, in the case of BH formation, POSYDON also al- 
lows us to take into account the detailed internal struc- 


ture of the star at the moment of collapse. As illustrated 
in the next section, we can then make an estimate of the 


spin of the CO taking into account the angular momen- 


tum profile of collapsing star. 


sin(0), where 0 


-falling material, j(r,9) = Qshen(r)r? 


in 


the polar angle, determines the properties of the ac- 
cretion flow. Disk formation occurs when the specific 


is 


angular momentum of the shell j(r, 0) exceeds the spe- 


cific angular momentum of the ISCO jīsco. This con- 


8.3.4. Birth spins of COs 


dition can be redefined as the polar angle at which disk 


formation occurs as 


We estimate the spin of the resulting BH following 
the collapse of the stellar profile as presented in Bavera 


(36) 


r 


( jisco 
Osheu ( 


arcsin 
The portion of the shell with 0 < 60a; will collapse 


Üaisk 


Appendix D). For convenience we summa- 


, 


(2021 
rize here the key assumption of this procedure. The final 


et al. 


mass and spin of the BH resulting from the collapse is 


calculated by following the accretion history of Mrembar 


directly onto the BH on a dynamical timescale trans- 
ferring j(r,0) to the hole, while the portion of the shell 


soon after the direct collapse of the central core of the 
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Figure 31. The CO mass Mco1, and spin aco: for the q = 0.7 grid slice of HMS-HMS MESA simulations in the initial primary 


mass-orbital period plane as predicted by the Patton & Sukhbold (2020) N20 engine. We assume that both the neutrino mass 
loss up to 0.5 Mo and the ejected mass during core collapse carries away the corresponding angular momentum. Spinning BHs 


are formed in binary systems avoiding mass transfer or undergoing stable mass transfer during contact phase or Case A mass 


transfer, see Figure 9. 


sionless spin parameter of the BH is updated after each 
shell is accreted onto the BH with the following relation 


with 0 > @gisk will form a disk and transfer only j1sco to 
the BH. The disk will be accreted on a viscous timescale 


which is assumed to be much smaller than the fallback 


(38) 


CJBH 


2 9? 
BH 


G 


a 


& Ramirez-Ruiz 


Batta 


timescale of the following shell ( 


2019). 


Therefore each collapsing shell contributes to the an- 


gular momentum of the BH by 


where Jpy is the angular momentum of the BH and 
Mpu its mass after accreting the directly infalling part 


the thin disk. 


The presented treatment is applicable only to the case 
of BH formation. For simplicity and the lack of firm 
alternatives in this version of POSYDON, we assume a zero 


spin for NSs. 


, 


of the shell and, if formed 


Jairect + Jaisk 
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Jshell 
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In the right panel of Figure 31, we show the CO spin 
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lar profile collapse assuming the remnant baryonic mass 
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of the ISCO of the accreting BH. This means that the 


shell COS(Üaisi 
1—e = 1 — [1 — 2GMgg/(3c rico)]!? is radiated away 


(Bardeen 1970; 


Thorne 1974 


, 


a 


is Maisk 


8.3.5. SN kicks 


fraction of the disk will be radiated away. The dimen- 


resultant BH will have mass smaller than Moray as a 
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During a SN, the binary system experiences abrupt 
mass loss, away from the center of mass, affecting 
its orbital parameters (Blaauw 1961; Boersma 1961). 
Furthermore, asymmetric ejection of matter (Janka & 
Mueller 1994; Burrows & Hayes 1996; Janka 2013) or 
asymmetric emission of neutrinos (Bisnovatyi-Kogan 
1993; Socrates et al. 2005) can provide a momentum 
kick to the newly formed CO. Here we assume that the 
magnitudes of the asymmetric kicks (vy) are drawn from 
a Maxwellian distribution with dispersion o: 


2 v? ie 
Fo) =y 775 exP (- 5) ; (39) 


As our fiducial assumption we take occgn = 265 km s7! 


(Hobbs et al. 2005) and egosw = 20kms' ! (Gia- 
cobbo & Mapelli 2019) for CCSN and ECSN, respec- 
tively. However, these velocities are free parameters. 
POSYDON supports multiple kicks rescaling options, e.g., 
if the prescription used to calculate the remnant bary- 
onic mass assume mass loss, i.e., the fallback mass frac- 
tion fg < 1, the kick is then rescaled by 1 — fm (Fryer 
et al. 2012). Alternatively BH kicks are rescaled by a 
factor 1.4 Mo /Mpu (using the gravitational mass) while 
NS kicks are not rescaled (our default option). Finally 
a user can opt to either not rescale any kicks or turn off 
SN kicks altogether. 

These kicks can tilt the orbit of the binaries, add ec- 
centricity or disrupt it. We take into account all these 
orbital changes including orbital changes for eccentric 
binaries following the analytical calculations of Kalogera 
(1996); Wong et al. (2012). 

We assume the collapsing star to lie on the origin of 
the coordinate system moving in direction of positive y- 
axis. The companion lies on the negative x-axis and z- 
axis completes the right-handed coordinate system (cf. 
Kalogera 1996, Figure 1). The semi-major axis after 
the kick a; is computed given the instantaneous orbital 
separation r; = a;[1 — ei cos(E)] pre-SN, where E is the 
eccentric anomaly, as 


vz + v? pue) (40) 


(2 
“=n GME, 


where M$ is the binary total stellar mass after the core 
collapse, v, is the pre-SN velocity of the collapsing star 
relative to the companion directed along the positive 
y-axis and v? the y-axis component of the kick. The 
eccentricity after the kick is then 


«-hi 


1/2 


(vk)? + [sin(v) (ve + vg) — cos(y)uk]? 2 


i i , 
GM iota 


where Mi, is the binary total stellar mass before the 
core collapse, w is the polar angle of the position vector 
of the collapsed star with respect to its pre-SN orbital 
velocity in the companion's reference frame and 


GM = e? )ai 


UrTi 


sin(v) = (42) 
In the above equations, in the case the CO receive no 
natal kick v = Okms-! but the star loses some mass, 
the orbit is still readjusted to conserve Kepler's third 
law. 

We consider a binary is disrupted if it does not sat- 
isfy the condition that demands the post-SN orbit passes 
through the pre-SN position (Willems et al. 2005; Flan- 
nery & van den Heuvel 1975), 


f. 
l-ep<—<lt+e, (43) 

af 
or if it is outside the limits of the amount of orbital 
contraction or expansion that can take place for a given 


amount of mass loss and a given magnitude of the kick 
velocity (Kalogera & Lorimer 2000; Willems et al. 2005) 


Mi Uk ? T; Mi Uk 2 
Qg——tet E A E E tot 1). 
Min C E ) Ag Mi Ur 
(44) 


Finally, we also verify that e; does not exceed 1 or that 
the argument of the square root in Eq. (41) does not be- 
come negative, if this is the case the binary is considered 
to be disrupted. 


9. HOW POSYDON EVOLVES AN INDIVIDUAL 
BINARY SYSTEM 


'To evolve a single binary within POSYDON, we use a hi- 
erarchy of classes. Every binary system is represented as 
a BinaryStar class containing two SingleStar classes, 
each with attributes that define their current state. To 
evolve the binary through each step, we have imple- 
mented a Pythonic flow, which takes the combination 
of a binary's state and event, and each stars’ state to 
direct a particular binary to its next step. We have a 
complete flow set as default, which can self-consistently 
track the evolution of binaries from their ZAMS state 
comprised of two H-stars, through all parts of the evo- 
lutionary tree shown in Figure 1. 

All steps in POSYDON are Python classes that update a 
binary via the user-defined call method. Steps are im- 
plemented based on on-the-fly calculations for evolution- 
ary phases such as the CE or the core collapse, or based 
on pre-calculated grids of detailed binary-star models, 
which is a novel component of POSYDON. In the latter 
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case, to estimate the evolution of systems for which no 
detailed model exists for the exact initial binary proper- 
ties, we use initial-final classification and interpolation 
algorithms, trained on the grids of detailed binary-star 
models (Section 7) or alternatively a nearest neighbor 
matching scheme. As the binary and its component 
stars evolve through these steps, the BinaryStar and 
SingleStar characteristics are appended to the objects, 
so every binary maintains a historical record of its evo- 
lution. 

If any system enters a phase that does not require fur- 
ther evolution, we use an end step to halt the evolution. 
This is used for binaries that either merge, disrupt, or 
reach the maximum physical time. For binaries whose 
evolution time ends in the middle of a step based on 
pre-calculated grid of models (e.g., if our star formation 
history randomly generates an initialization time for a 
binary within a few megayear of the end run time), we 
can no longer use our pre-trained classification and in- 
terpolation algorithms as these only apply to the end 
state of the binary. In these cases, we instead use the 
system's nearest-neighbor pre-computed track directly 
(interpolating full binary tracks is non-trivial and is be- 
ing investigated for future versions of POSYDON). This is 
the default behavior of all MESA grid steps, with various 
classification and interpolation algorithms ready to use. 

Figure 32 depicts the complete evolution of one par- 
ticular binary from ZAMS to the formation of a BBH 
system that merges within a Hubble time. Each vertical, 
colored band indicates an extended stage of a binary's 
evolution, whereas the two core-collapse events and the 
CE phase are essentially instantaneous processes, occur- 
ring in between the other, extended phases. 

The initial masses of the system are M, = 97.07 Mo, 
Mz = 39.86 Moin a Py = 39.45 days circular orbit. 
The first part of the evolution of the system is based on 
the HMS-HMS binary grid (Section 5.5), and its sub- 
sequent evolution is followed either through the nearest 
neighbor interpolation (orange lines) or through initial- 
final interpolation (black dots; using our classification 
and interpolation methods) for the same initial configu- 
ration. In both cases we adopt the same SN kick after 
the two core-collapse episodes in order to compare them 
as close as possible. 

'he system does not experiences mass transfer be- 
fore the first supernova, however the primary experi- 
ence strong stellar wind mass loss during the Wolf-Rayet 
phase which widens the orbit . Shortly after the primary 
collapses into a ~ 14.96 Mo BH with a low (aspin),- 
The subsequent detached evolution (Section 8.1) of the 
mildly eccentric system (due to the BH natal kick), after 
matching the companion of the BH to a single star grid, 


leads to a small increase of the period predominantly due 
to winds. Eventually the secondary star fills its Roche 
lobe at periastron, where we assume that the system 
circularizes. Mass transfer onto the BH is interpolated 
through the binary grid of COs with H-star companions 
(Section 5.6) and lasts for a few thousand years, becom- 
ing unstable and leading to a CE episode (Section 8.2). 
'The system survives the process, forming a tight binary 
comprised of a BH with a stripped He star on a ~ 0.2 day 
orbit. Tidal forces become important in this tight orbit, 
spinning up the He star (Section 5.7), which eventu- 
ally also forms a mildly spinning (dspin,2 ~ 0.48) BH of 
1.90 Mo. The two BHs merge after 183 Myr from birth, 
due to gravitational wave radiation. 

We have specifically chosen a binary where the dif- 
ferences between the nearest neighbor and initial-final 
interpolation schemes are relatively small, so that we 
can accurately display how the binary evolves through 
each step. Differences between the two evolution options 
for binaries in general are significantly larger, and the 
initial-final interpolation scheme is our default choice. 

'The binary shown in Figure 32 was evolved using our 
default configuration, although we have purpose-built 
POSYDON to be modular. Throughout the previous sec- 
tions we have described possible changes to physical pre- 
scriptions that a user can make. However, a user can also 
easily supplement their own functions for specific steps, 
or even define the entire binary flow. 

Finally, for debugging we keep track of any errors or 
warnings raised throughout a binary’s evolution. This 
allows us to isolate the problematic step for a binary, or a 
stellar and binary state-event combination that our flow 
structure cannot handle. This error tracking can be es- 
pecially useful for user-defined steps and flow structures. 


10. HOW POSYDON EVOLVES A BINARY 
POPULATION 


Generating a model binary population for compari- 
son to observations requires two separate steps: initial- 
izing individual binaries and then evolving those bina- 
ries, which we describe in Section 10.1 and Section 10.2 
below. Finally, in Section 10.3 we describe a sample 
binary population evolved with POSYDON. 


10.1. Initialization 


'The primary function of POSYDON is to produce syn- 
thetic populations of binaries which requires evolving a 
random distribution of binaries from ZAMS. We gener- 
ate an initial population by sampling binary parameters 
from standard distributions. 

Component Masses: For the primary component 
mass of a binary, we implement initial mass func- 
tions (IMF) from Salpeter (1955); Kroupa et al. (1993); 
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Figure 32. Evolution over time, from ZAMS to binary BH formation, of one example binary system with initial properties 
Mi = 97.07 Mo, M» = 39.86 Mo in a Pa» = 39.45 days. The top two rows of panels show the evolution of the binary's Por» and e 
for our nearest neighbor matching scheme (orange lines) and our initial-final interpolation method (black, circle markers). In the 
bottom three rows, we show the primary star's (1; solid lines) and secondary star's (2; dashed-dotted lines) properties using the 
nearest neighbor interpolation method and compare them against the same binary using the initial-final interpolation method 
(circle markers for primary, star markers for secondary). The binary's evolution is followed across its different evolutionary steps 
(note that the time scale varies for each step), through both MESA grids (Section 5) and on-the-fly calculations (Section 8). 


Kroupa (2001). By default we use the Kroupa (2001) 
IMF in the range Mı € [7,150] Mo with a = 2.3. Then 
the secondary component mass is given by drawing the 
mass ratio q = M2/M, from a flat distribution [0, 1] with 
M» € [0.35, 150] Mo. 

Orbital Parameters: For a binary population's ini- 
tial orbital period, Sana et al. (2013) describe a power- 
law in log-space, while Opik's Law describes a log-flat 
distribution in orbital separation rather than orbital pe- 
riod space (Abt 1983). 

In POSYDON we allow for both models to be adopted, 
with our default being Sana et al. (2013). The mini- 
mum orbital period is set by systems that undergo RLO 
at ZAMS, and a default maximum orbital period to be 


6000days. The orbital period distribution from Sana 
et al. (2013) is undefined for Por» < 1day; therefore, we 
extend this distribution so that it is uniform in logy Porb 
space down to a Py of 0.35 days. We further set the 
maximum P, to 10^? days. Although binaries may 
form at wider separations, we choose this upper limit to 
account for all interacting binary models (e.g., Figure 9). 

Since we only model circular binaries in our detailed 
binary-star models with the MESA code, the binary ec- 
centricity e at ZAMS is set to zero (see Section 4 and 
Section 8.1). Since we plan to generalize this assumption 
in future work, we also add the option to generate e from 
a thermal distribution (Duquennoy & Mayor 1991). 
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Figure 33. Example of the double compact object populations produced by POSYDON. Quantities shown are the final distribu- 
tions, as the binary populations appear today. We separately indicate the parameters for NS-NS, BH-NS, and BH-BH systems. 
Results from a population synthesis using linear interpolation is indicated with solid lines, while results using 
the first nearest neighbor approach with dashed lines. Comparison with COSMIC (Fig. 3 from Breivik et al. 2020) shows 
morphological similarities, but important quantitative differences between the same binary CO populations. 


Star-formation History: To account for the star- 
formation history (SFH) of a binary population, we as- 
sign to each binary in the population the same maxi- 
mum age. Then different star-formation rates (SFRs) 
can be modeled by modifying the distribution of birth 
times. In POSYDON we offer two options for the SFR: a 
burst of star formation or a constant SFR. For the burst 
model, all binaries have identical birth times some num- 
ber of years prior to the end of the simulation. For a 
constant SFR, we randomly generate birth times from 
a uniform distribution within a user-defined range. The 
average metallicity of the Galaxy and the greater Uni- 
verse evolve over time, but this is something that we 
cannot currently model accurately, as our grids of mod- 
els presented in v1.0 of POSYDON are only calculated for 
stars at solar metallicity. 


10.2. Evolution 


Evolving a population of ZAMS binaries initialized us- 
ing the distributions described in Section 10.1 requires 
implementing the procedure outlined in Section 9 for 
each binary. In POSYDON we have created an overarching 
BinaryPopulation class which is a container for a list 
of individual BinaryStar instances. Each BinaryStar 
instance is then iteratively evolved until the entire pop- 
ulation has been processed. 

The BinaryPopulation class contains a number of ad- 
ditional capabilities to efficiently and easily evolve pop- 
ulations of binaries. First, evolved populations are au- 
tomatically saved in an efficient hdf5 file format with 
two datasets: one that contains each binary as a single 
line providing both its initial and final states, and a sec- 
ond dataset that contains the entire evolutionary history 
of each binary. Second, populations can be evolved ei- 
ther serially or in parallel, so that large (> 10° binaries) 


populations can be run quickly on a high-performance 
computing cluster. Third, we have implemented rou- 
tines that catch and keep warnings and errors from each 
binary, so that the code does not crash when a single 
binary fails to complete its evolution due to a bug. We 
have found this to be useful for identifying and resolving 
coding bugs when implementing new physical prescrip- 
tions. Last, we have implemented various routines that 
allow a user to easily select only certain types of binaries 
(e.g., only BBHs or only double COs). 

We envision that a typical user will interact primar- 
ily with the BinaryPopulation class and its associ- 
ated SimulationProperties class, which when com- 
bined provide the interface for customizing a particular 
user’s binary population synthesis needs. 


10.3. Example Population 


'To demonstrate the results of a binary population syn- 
thesis run with POSYDON, we construct a basic popula- 
tion of 10° binaries generated with the default initial 
conditions described in Section 10.1. We choose a con- 
stant star formation history over the past 10 Gyr. The 
simulation on our high performance computing cluster, 
Trident, takes approximately 2 hours of wall time using 
5 nodes, each with 20 cores, i.e., less than 1s of CPU 
time per binary. 

From the resulting binaries, we select those that evolve 
into bound NS-NS$S, NS-BH, or BH-BH systems. The 
present-day properties of these binaries are provided in 
Figure 33; this can be compared with the results of other 
studies, e.g., Figure 3 from Breivik et al. (2020) with 
similar initial conditions at solar metallicity. Overall we 
find good agreement reaffirming that the code is produc- 
ing reasonable results. In the left two panels, we show 
distributions of the component compact object masses. 
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Our NS masses are in very close agreement with those of 
Breivik et al. (2020), although the BHs we produce ex- 
tend to larger masses. The reason is that current stellar 
models with stellar structure and evolution parameters 
calibrated to the latest observations (e.g., overshoorting, 
etc) produce cores more massive compared to those in 
the late-nineties models used in rapid BPS codes. The 
third panel from the left shows the semi-latus rectum of 
the population. Again, the distribution results are qual- 
itatively similar, with a preference for NS-NS systems 
at smaller a(1 — e”), and NN-BH and BH-BH systems 
having larger a(1 — e?). Finally, the rightmost panel of 
Figure 33 shows an increasing formation timescale (from 
ZAMS to the second SN) as we move from BH-BH to 
NS-BH to NS-NS systems. This is expected since BHs 
tend to form from more massive systems that complete 
their evolution more quickly. Our binary CO popula- 
tions appear morphologically similar to those produced 
by COSMIC and described in Breivik et al. (2020). We 
will undertake detailed descriptions of the specific bi- 
nary populations of interest to separate, science studies. 


11. SUMMARY AND FUTURE WORK 


Here we present POSYDON, a new, next-generation com- 
putational tool for general population synthesis of single 
and binary stars. POSYDON incorporates full stellar struc- 
ture and evolution sequences for interacting binaries, us- 
ing the MESA code. Compared to other existing, binary 
population synthesis code there are significant advances: 
(i) binary evolution is treated self-consistently without 
analytical fits of single-star evolutionary tracks and the 
need for simplified or artificial recipes to emulate the 
behavior of stars in interacting binaries; (ii) initial-final 
classification and interpolation methods trained on the 
pre-calculated grids of binary evolution models, allow- 
ing general synthetic simulations of binary populations. 
'The code base along with the existing evolutionary-track 
grids are publicly available through the POSYDON collab- 
oration's web portal (https://posydon.org) along with 
full documentation and tutorials for how to use the 
code. An advanced query system is also available for 
users to be able to mine the grids of single- and binary- 
star evolutionary tracks and download relevant data us- 
ing pre-programmed and customized queries (e.g., Teng 
et al. 2021a,b). Finally, we provide a user-friendly web- 
application that allows a user to perform small-scale 
simulations with POSYDON online, without the need of 
code installation and configuration. 

Compared to current rapid population synthesis codes 
POSYDON has a smaller set of free parameters, for many 
of which there are already multiple options for the user 
to choose from. The code structure is modular and an 


advanced user is able to implement their own choices 
of evolutionary parameters; from as simple as changing 
the initial properties of the binaries to as complex as in- 
corporating their own custom-made evolutionary-track 
grids. In this first instrument POSYDON paper we describe 
in detail the first version of the code, but technical and 
astrophysical advancements are ongoing and improved 
code versions will be released in the near and long-term 
future. 

Our focus on the technical front is on classification and 
interpolation methods. Our current process of first clas- 
sifying the grids and then performing interpolation can 
lead to errors propagating throughout the pipeline. In- 
stead these two could be combined into a joint treatment 
to reduce errors, making additional use of covariances 
between the different grid types (Singh et al. 2016). We 
will also explore adopting kernelized-interpolation ap- 
proaches (i.e., Wilson & Nickisch 2015; Gardner et al. 
2018; Narayan et al. 2021). To illustrate the use of such 
an approach, it is observed in Figure 22 that the con- 
fusion matrix for our grid of He-rich stars with com- 
pact object companions has a lower accuracy for un- 
stable mass transfer. This is perhaps due to the non- 
linearly separable decision boundary, as visible in Fig- 
ure 21. Adopting a kernelized-interpolation technique 
has the potential to increase the class-specific interpo- 
lation accuracy. In particular, kernel-based approaches 
such as support vector machines and Gaussian processes 
may prove fruitful. In principle, we may also use a neural 
networks. Also, improvements can be made by adopting 
non-Euclidean metrics when defining our distance func- 
tions in our k-nearest neighbors classifiers in Section 7.2. 

Apart from methods exploration for the existing clas- 
sification and interpolation processes, we will focus on 
the next step of interpolations critical for astrophysical 
studies: interpolation between whole evolutionary tracks 
along time. This is a challenging problem that is of par- 
ticular interest for any study that requires tracking bi- 
nary properties as a function of their age (e.g., X-ray 
binary luminosities). In parallel, we are already work- 
ing on increasing the computational efficiency of build- 
ing the pre-computed grids necessary for future POSYDON 
versions. Specifically, for any future grid development, 
we will take advantage of a new active-learning method 
developed by our team (Rocha et al. 2022) that allows 
us to achieve the same classification and interpolation 
accuracies with a significantly smaller MESA tracks, by 
dynamically placing them at class boundaries in the 
parameter space. Such dynamic placement, informed 
by active learning, leads to grids with non-regular, but 
smart placement of evolutionary tracks. This work be- 
comes more critical as we expand to more metallicities 
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and add eccentricity as another dimension in order to 
keep the POSYDON package sizes appropriate for down- 
loads. 

The first version of POSYDON is fully functional as an 
astrophysical tool in the sense that it can be used for 
complete simulations of binary populations from ZAMS 
to either formation of a binary with two compact objects 
or binary distractions (stellar merger or binary disrup- 
tion), but it is still limited in two ways: our current 
pre-computed grids are for primaries massive enough to 
likely form a NS or a BH and are calculated at solar 
metallicity only. Our next version will expand to a grid 
of metallicities appropriate for populations across the 
Universe instead of just the Milky Way. Challenges re- 
lated to the convergence of single- and binary-star evo- 
lutionary models often depend on metallicity, e.g. be- 
cause the different opacities will results in massive stars 
reaching the Eddinghton limit within their interiors at 
different mass ranges and evolutionary phases. Typi- 
cally, however, such challenges become less severe with 
decreasing metallicity. Interpolation across metallicities 
will be a follow-up step as well, once grids at a sufficient 
number of metallicities have been computed. All the 
interpolation and classification methods we have devel- 
oped scale naturally to additional dimensions in the ini- 
tial conditions parameter space. Further we will expand 
our grids to low-mass primaries so populations with 
WDs can also be modeled. Expansion of the grids will 
inevitably increase data set sizes proportionally. The de- 
velopment of interpolation methods for not only the final 
properties of each grid but also the entire evolutionary 
tracks, which is one of our primary future objectives, 
will allow future POSYDON releases to come only with 
with the pre-trained interpolation objects. The latter 


are expected to have a significantly smaller data foot- 
print compared to the downsampled grids we currently 
ship with POSYDON. 

We will also continue to improve the physics treat- 
ment of binary evolution. Specifically we are already 
working on three improvements: (i) in this first version 
our treatment of the CE phase, once dynamical insta- 
bility is recognized (taking into account the full stellar 
structure of the RLO star), is similar to what is done in 
pBPS codes, apart from the self-consistent calculation of 
the CE's binding energy. However, since we model bina- 
ries with MESA, we are able to treat the phase in a more 
physical way, either by following the CE inspiral self- 
consistently using one-dimensional hydrodynamic sim- 
ulations (e.g., Fragos et al. 2019) or by following the 
long-term response of the RLO star to losing its enve- 
lope on a very high rate and have its Roche lobe shrink- 
ing rapidly (e.g., Marchant et al. 2021; Gallegos-Garcia 
et al. 2021). One of the objectives of the next version 
of POSYDON will be to improve the physics of our CE 
treatment. (ii) like all binary population modeling to 
date, we assume that binaries circularize instantly upon 
RLO and for this reason we also assume that all ZAMS 
binaries are circular. However, the physics of secular bi- 
nary evolution through mass transfer in eccentric orbits 
has been fully developed recently (Sepinsky et al. 2007, 
2009, 2010; Dosopoulou & Kalogera 2016a,b; Hamers 
et al. 2021) and will be implemented in future POSYDON 
versions. (iii) Most recently more physical models for 
magnetic braking have been developed and calibrated 
against single-star rotational-velocity data (e.g., Van & 
Ivanova 2019, 2021; Gossage et al. 2021) and we will 
use them to update the options for magnetic braking 
evolution in binaries. 


APPENDIX 


'Table 5. The variables in history tables in our grids, taken from the MESA output. The star.age is not included since we provide 
the age variable in the binary history (cf. Table 6) which refers to both the systems, and its components. Since these quantities are 
undefined for compact objects, history tables are not provided for NSs and BHs; however, the evolution of the mass is given in the 
binary history tables. 


Name Description Unit 
he.core.mass Helium core mass. Mo 
c_core_mass Carbon core mass. Mo 
o_core_mass Oxygen core mass. Mo 
he_core_radius* Helium core radius. Ro 
c-core.radius * Carbon core radius. Ro 


Table 5 continued 
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Table 5 (continued) 


Name Description Unit 
o.core.radius * Oxygen core radius. Ro 
center hl Center ! H mass fraction. 
center he4 Center ^He mass fraction. 
center. c12 * Center !2C mass fraction. 
center. n14 * Center 14N mass fraction. 
center. 016 * Center 160 mass fraction. 
surface. h1* Surface !H mass fraction. 
surface he4 * Surface ^He mass fraction. 
surface. c12* Surface 12C mass fraction. 
surface.n14* Surface !^N mass fraction. 
surface.o16 * Surface 160 mass fraction. 
c12.c12* Decimal logarithm of the burning power from the !?C + !1?C reaction [Lo] 
center_gamma* Plasma coupling parameter, ratio of the Coulomb to thermal energy. 
avg.c.in.c.core Average !1?C abundance at carbon core. 
surf avg omega * Average surface angular velocity. yr 
surf.avg-omega.div.omega.crit * Ratio of the average and critical surface angular velocity. 
log. LH * Decimal logarithm of the hydrogen burning power. Lo 
log-LHe* Decimal logarithm of the helium burning power. Lo 
log.LZ* Decimal logarithm of the total burning power excluding LH and LHe and [Lo 

photodisintegration. 
log.Lnuc * Decimal logarithm of the total nuclear burning power. Lo 
log-Teff Decimal logarithm of the effective temperature. K] 
log.L Decimal logarithm of the luminosity. Lo 
log_R Decimal logarithm of the radius. Ro 
total moment. of inertia 'Total momentum of inertia. gcm? 
spin. parameter * Dimensionless stellar spin parameter. 
log.total.angular. momentum Decimal logarithm of the total angular momentum. gcm? s-1] 
conv.env.top.mass * Mass coordinate of the top boundary of the outermost convective region. Mo 
conv.env. bot. mass * Mass coordinate of the bottom boundary of the outermost convective region. Mo 
conv.env. top. radius * Radial coordinate of the top boundary of the outermost convective region. Ro 
conv.env.bot. radius * Radial coordinate of the bottom boundary of the outermost convective region. Ro 
conv_env_turnover_time_g* Global convective turnover time. yr 
conv.env.turnover.time.l b * Local convective turnover time half of a scale height above the outermost convective yr 
zone bottom boundary. 
conv.env.turnover.time.l t * Local turnover time one scale height above the outermost convective zone bottom yr 
boundary. 
envelope. binding. energy * Binding energy of the envelope. erg 
mass. conv. reg. fortides * Mass of the most important convective region for equilibrium tides, as defined in Mo 
Eq. (7). 
thickness. conv.reg. fortides * Thickness of the most important convective region for equilibrium tides, as defined Ro 
in Eq. (7). 
radius.conv.reg.fortides * Radial coordinate of the most important convective region for equilibrium tides, as Ro 


lambda. CE. 1cent * 
lambda. CE. 10cent * 


lambda. CE. 30cent * 


defined in Eq. (7). 


Common-envelope parameter of the envelope binding energy for core-envelope bound- 
ary where hydrogen mass fraction becomes lower than 196. 


Common-envelope parameter of the envelope binding energy for core-envelope bound- 
ary where hydrogen mass fraction becomes lower than 10%. 


Common-envelope parameter of the envelope binding energy for core-envelope bound- 
ary where hydrogen mass fraction becomes lower than 30%. 


'Table 5 continued 
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Name Description Unit 
co_core_mass Carbon-oxygen core mass. Mo 
co_core_radius* Carbon-oxygen core radius. Ro 
lambda. CE. pure. He. star. 10cent * Common-envelope parameter of the He-rich envelope binding energy for core-envelope 

boundary where the sum of hydrogen and helium mass fraction becomes lower than 
10%. 
log_L_div_Ledd* Decimal logarithm of the ratio of the luminosity and Eddington luminosity. 
*Property not accounted for when downsampling the grids (Section 6.4). 
Table 6. The variables in binary history tables, taken from the MESA output. 
Name Description Unit 
model. number * 'The model number of the final state 
age" Binary age yr 
star.1 mass Mass of the first star Mo 
star_2_mass Mass of the second star Mo 
period_days Orbital period in days d 
binary_separation Binary separation Ro 
g-system. mdot. 1 Decimal logarithm of rate of mass loss from the system from around the first star [Mo yr"! 
due to inefficient mass transfer 
g-system. mdot. 2 Decimal logarithm of rate of mass loss from the system from around the second star [Mo yr! 
due to inefficient mass transfer 
g_wind_mdot_1* Decimal logarithm of rate of mass loss of the first star due to wind Mo yr 
g_wind_mdot_2* Decimal logarithm of rate of mass loss of the second star due to wind Mo yr 
g_mstar_dot_1* Decimal logarithm of rate of mass loss of the first star Moyr7t 
g_mstar_dot_2* Decimal logarithm of rate of mass loss of the second star Mo yr- 
g.mtransfer. rate Decimal logarithm of mass-transfer rate Mo yr 
xfer_fraction* Mass-transfer fraction 
rl_relative_overflow_1* Roche lobe overflow of the first star in units of donor Roche lobe radii 
rLrelative_overflow_2* Roche lobe overflow of the second star in units of donor Roche lobe radii 
trap_radius* Trapping radius Ro 
acc_radius* Radius of the compact object cm 
t.sync.rad.1* Tidal synchronization time-scale of the first star for stars with radiative envelopes S 
t_sync_conv_1* Tidal synchronization time-scale of the first star for stars with convective envelopes s 
t_sync_rad_2* Tidal synchronization time-scale of the second star for stars with radiative envelopes s 
t_sync_conv_2* Tidal synchronization time-scale of the second star for stars with convective envelopes s 
*Property not accounted for when downsampling the grids (Section 6.4). 
Table 7. Quantities of final profiles of the stars, taken from the MESA output. 
Name Description Unit 

radius Radius at the outer boundary of the zone. [Ro] 

mass* Mass coordinate of the outer boundary of the zone. [Mo] 

logRho Decimal logarithm of the density at the center of the zone. [g cm 3] 
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Table 7 (continued) 


Name 


Description 


Unit 


omega 
energy j 
x_mass_fraction_H* 


y-mass_fraction_He i 


z mass fraction metals * 


neutral fraction H * 
neutral fraction. He* 


avg. charge He * 


Angular velocity. 


Specific internal energy. 


[rads7 1] 


[erg g- 1] 


Mass fraction of all isotopes with atomic number 1. 


Mass fraction of all isotopes with atomic number 2. 


Mass 


fraction of all 


elements except for those in x mass fraction H and 


y-mass.fraction. He 
Fraction of neutral hydrogen (HI) of all the !H. 
Fraction of neutral helium (Hel) of all the 4He. 


Average charge of all the ^He isotopes. 


electron charge [e] 


* Property not accounted for when downsampling the grids (Section 6.4). 


Table 8. Post-processed variables referring to the final state of each binary in our grids. Quantities referring to single- 


star quantities (i.e., all except termination flags and interpolation class) are prefixed with S1. and S2. to distinguish the 


corresponding star in the grid (e.g., S1. surface other or 82 direct mass). These variables, along with the last values of the 


single and binary history variables (cf. Table 5 and Table 6; e.g., S1 log L), comprise the final values tables stored in the 


grids. 

Name Description Unit 
termination. flag 1 Termination reason from MESA output, or reach cluster timelimit. 
termination_flag_2 RLO state (indicating which star is the donor), or contact during MS in case 

of stellar merger. 
termination. flag. 3 State of primary star. 
termination flag 4 State of secondary star. 
interpolation. class Classification based on termination flags 1 and 2, indicating broad groups 
based on mass transfer. 
surface. other Surface abundance fraction of elements excluding !H, He, !2C, 14N, 160. 
center_other Central abundance fraction of elements excluding 'H, ^He, !2C, 14N, 160, 
direct. state CO state for the direct collapse prescription of the star. 
direct. SN type SN type for the direct collapse prescription of the star. 
direct. f fb Fallback mass fraction for the direct collapse prescription of the star. 
direct. mass CO mass for the direct collapse prescription of the star. Mo 
direct.spin CO spin for the direct collapse prescription of the star. 
Fryer-4-12-rapid state CO state for the Fryer et al. (2012) rapid prescription of the star. 
Fryer4-12-rapid. SN type SN type for the Fryer et al. (2012) rapid prescription of the star. 
Fryer4-12-rapid f fb Fallback mass fraction for the Fryer et al. (2012) rapid prescription of the 
star. 
Fryer--12-rapid. mass CO mass for the Fryer et al. (2012) rapid prescription of the star. Mo 
Fryer+12-rapid_spin CO spin for the Fryer et al. (2012) rapid prescription of the star. 
Fryer4-12-delayed.state CO state for the Fryer et al. (2012) delayed prescription of the star. 
Fryer4-12-delayed.SN type SN type for the Fryer et al. (2012) delayed prescription of the star. 
Fryer4-12-delayed. f fb Fallback mass fraction for the Fryer et al. (2012) delayed prescription of the 
star. 
Fryer4-12-delayed mass CO mass for the Fryer et al. (2012) delayed prescription of the star. Mo 


Fryer4 


-12-delayed.spin 


Sukhbold-F-16-engineN20.state 
Sukhbold--16-engineN20.SN. type 


CO spin for the Fryer et al. (2012) delayed prescription of the star. 
CO state for the Sukhbold et al. (2016) N20 engine prescription of the star. 
SN type for the Sukhbold et al. (2016) N20 engine prescription of the star. 


'Table 8 continued 
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Table 8 (continued) 


Name Description Unit 
Sukhbold+16-engineN20_f fb Fallback mass fraction for the Sukhbold et al. (2016) N20 engine prescription 
of the star. 
Sukhbold+16-engineN20_mass CO mass for the Sukhbold et al. (2016) N20 engine prescription of the star. Mo 
Sukhbold--16-engineN20.spin CO spin for the Sukhbold et al. (2016) N20 engine prescription of the star. 
Patton&Sukhbold20-engineN20. state CO state for the Patton & Sukhbold (2020) N20 engine prescription of the 
star. 
Patton&Sukhbold20-engineN20. SN type SN type for the Patton & Sukhbold (2020) N20 engine prescription of the 
star. 
Patton&Sukhbold20-engineN20.f fb Fallback mass fraction for the Patton & Sukhbold (2020) N20 engine pre- 
scription of the star. 
Patton&Sukhbold20-engineN20. mass CO mass for the Patton & Sukhbold (2020) N20 engine prescription of the Mo 
star. 
Patton&Sukhbold20-engineN20. spin CO spin for the Patton & Sukhbold (2020) N20 engine prescription of the 
star. 
avg.c.in.c.core.at. He. depletion Average carbon 12 abundance at carbon core at the state of He depletion of 
the star. 
co.core.mass.at. He.depletion Carbon-oxygen core mass at the state of He depletion of the star. 
m.core. CE. 1cent Mass of the hydrogen-deficient (i.e. helium) core, with the core-envelope Mo 
boundary defined as the outermost layer where the hydrogen mass fraction 
drops below 196. 
m. core. CE. 10cent As m_core_CE_lcent, but for hydrogen mass fraction of 10%. Mo 
m_core_CE_10cent As m_core_CE_lcent, but for hydrogen mass fraction of 30%. Mo 
m_core_CE_pure_He-star_10cent Mass of the hydrogen- and helium-deficient (i.e., carbon-oxygen) core, with Mo 
the core-envelope boundary defined as the outermost layer where the sum of 
hydrogen and helium mass fraction drops below 1096. 
r core CE. 1cent Radial coordinate of the core as defined in m. core CE. 1cent. Ro 
r_core_CE_10cent Radial coordinate of the core as defined in m_core-CE_10cent. Ro 
r core CE 30cent Radial coordinate of the core as defined in m. core. CE. 30cent. Ro 
r_core_CE_pure_He-_star_10cent Radial coordinate of the core as defined in m_core.CE_pure_He-_star_10cent. Ro 
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